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Preface 


This  reporr  presents  an  evaluation  of  mass  Integral  scaling  (MIS) 
as  applied  to  the  atmospheric  radiation  transport  problem.  MIS  is  the 
fundamental  approximation  used  In  scaling  codes  like  SMAUG  and  ATR  to 
predict  fnee-field  radiation  environments  from  atmospheric  nuclear 
detonations.  Such  predictions  are  required  in  and  are  an  integral  part 
of  analyses  designed  to  assess  the  nuclear  survivability/vulnerabi  lUy  of 
military  equipment,  systems,  and  personnel. 

The  results  of  the  evaluation  are  summarized  in  a series  of  iso- 
error contour  graphs  which  clearly  define  the  MIS  associated  error  over 
various  spatial  domains  in  the  1-25  km  altitude  region.  These  contours 
show  that,  while  MIS  is  an  acceptable  approximation  for  source  altitudes 
between  1 and  10  km,  errors  of  a factor  of  two  or  more  are  easily 
incurred  if  applied  at  higher  altitudes. 

The  purpose  of  the  research  reported  here  is  not  to  condeim  the  use 
of  the  popular,  quick-running  scaling  codes.  Instead,  iry  intent  was  to 
qualify  the  results  of  these  codes.  In  fact,  I believe  that  the  codes 
are  even  more  valuable  now,  because  for  the  first  time,  users  know  how 
bad  or  good  their  results  are.  Furthermore,  I believe  that  an  MIS 
correction  factor  algorithm  can  be  fairly  easily  constructed  using  the 
data  and  methodology  developed  in  the  course  of  this  study. 

I would  like  to  acknowledge  several  people  who  provided  me  with 
significant  assistance  during  the  course  of  this  effort.  First,  I want 
to  thank  nqy  dissertation  advisor.  Dr.  Charles  J.  Bridgman  (Air  Force 
Institute  of  Technology),  for  his  suggestions  and  the  encouragement 
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that  he  provided.  He  was  one  of  the  first  to  question  the  validity  of 
HIS  and  convinced  me  that  an  evaluation  of  the  technique  was  essential. 
Next,  1 want  tu  express  ny  appreciation  to  ny  friends  and  colleagues, 

John  Burgio  and  Harry  Murphy  of  the  Air  Force  Weapons  Laboratory.  John 
processed  the  cross  sections  used  in  the  study  and,  during  our  many 
discussions,  provided  me  with  a much  broader  insight  to  the  problem. 

Harry  deserves  special  mention.  He  Is,  without  a doubt,  the  most  com- 
petent computer  programmer  that  J.  have  ever  met.  He  wrote  the  fitting 
programs  used  in  the  study  and  provided  me  with  several  routines  which 
were  incorporated  directly  into  the  AIRDIF  code.  The  structured,  modular 
nature  of  AIRDIF  is  attributable  to  the  patient  instruction  he  gave  me. 

I am  also  indebted  to  Mrs.  Kay  Fink  and  Ms.  Maria  Garcia  for  typing  the 
drafts  of  this  report,  to  Mrs.  Bette  Shannon  for  typing  the  final  manu- 
script, to  Mr.  Clarence  Larson  for  the  preparation  of  all  figures  and 
graphs,  and  to  ny  math  aide,  Keith  Merrifield,  for  the  many  hours  he 
spent  in  helping  me  polish  the  various  drafts  into  final  form.  Finally, 

I want  to  t iank  ny  family  and  especially  ny  loving  wife,  Vicki,  for  the 
patience  and  understanding  they  showed  me.  Their  sacrifices  over  the 
past  two  years  were  far  greater  than  mine. 

Raymond  A.  Shulstad 
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Abstract 


This  report  presents  the  results  of  a comprehensive  evaluation  of 
mass  integral  scaling  (MIS)  as  applied  to  the  atmospheric  radiation 
transport  problem.  MIS  *s  an  approximation  used  to  predict  the  free- 
field  radiation  environments  from  atmospheric  nuclear  detonations. 

Codes  like  SMAUG  and  ATR  use  macs  integral  scaling  of  infinite  homoge- 
neous air  transport  data  to  generate  environments  which  purport  to 
reflect  the  near  exponential  vertical  variation  in  density  of  air. 

The  validity  of  the  MIS  approximation  was  evaluated  by  determining 
K factors  defined  to  be  the  ratio  of  the  2-D  4ttR2  variable  density  air 
dose  to  the  1-D  4ttR2  MIS  homogeneous  air  dose.  K factors  were  determined 
for  neutron  and  secondary  gamma  doses  from  both  an  unclassified  fission 
and  a thermonuclear  source.  Calculations  were  made  for  source  altitudes 
between  5 and  20  km,  for  receiver  altitudes  between  1 and  25  km,  and  for 
mass  ranges  out  to  250  g/cm*.  The  results  are  presented  in  the  form  of 
a series  of  K factor  versus  mass  range  curves  for  a particular  source 
spectrum  and  source/receiver  altitude.  These  curves  may  be  used  to 
extract  multiplicative  MIS  correction  factors.  A series  of  MIS  iso- 
error contours,  which  clearly  define  the  error  bounds  over  various 
spatial  domains,  is  also  presented. 

The  unique  and  distinguishing  feature  of  this  research  is  that  the 
evaluation  was  conducted  with  a special  form  of  diffusion  theory  which 
is  also  described  in  this  report.  By  comparing  the  diffusion  derived 
K factors  to  some  recent  Monte  Carlo  calculations,  it  is  shown  that  the 
evaluation  is  not  dependent  on  diffusion  theory.  This  is  significant 
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because  the  full-scale  parametric  evaluation  of  MIS  was  conducted  with 
diffusion  theory  at  a fraction  of  the  cost  that  would  have  been  required 
by  a comparable  Monte  Carlo  evaluation. 

The  results  of  this  evaluation  indicate  that  MIS  is  acceptable 
(errors  less  than  20%)  out  to  mass  ranges  of  220  g/cm2  for  source 
altitudes  between  1 and  10  km.  At  higher  altitudes,  however,  errors 
of  a factcr  of  two  or  more  are  easily  incurred. 
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AN  EVALUATION  OF  MASS  INTEGRAL  SCALING 
AS  APPLIED  TO  THE 

, ATMOSPHERIC  RADIATION  TRANSPORT  PROBLEM 

I.  Introduction 

Background 

The  prediction  of  f ree-fie Id  radiat  ion  environments  from  atmospneric 
nuclear  detonations  is  an  integral  part  of  analyst-  designed  to  assess 
the  survivability/ vulnerability  (S/V)  of  military  equipment*  systems,  and 
personnel.  These  assessments  are  characterized  by  the  requirement  to 
evaluate  survivability  and  vulnerability  under  a multitude  of  scenarios. 
For  example*  the  analyst  must  consider  the  sensitivity  of  his  results  to 
such  parameters  as  weapon  yield,  source  spectrum,  and  burst/ receiver 
altitudes. 

The  physics  of  the  atmospheric  radiation  transport  problem  is 
described  by  the  well-known  Boltzmann  transport  equation.  Unfortunately, 
because  of  the  unique  problems  created  by  the  exponential -like  density 
variation  in  air*  (Figure  1),  classical  numerical  solutions  of  the 
Boltzmann  equation  for  this  problem  have  not  been  attempted.  Monte 
Carlo  simulation  is  the  only  method  that  has  been  used  to  compute  the 
free-field  radiation  environments  in  two-dimensional  (2-D)  variable 
density  air.  Table  I summarizes  the  limited  number  of  2-D  air  Monte 
Carlo  calculations  which  have  been  published. 

♦Variable  density  air  In  this  study  is  defined  to  be  the  1962  U.S. 
Standard  Atmosphere  (Pef.  1). 


1 


TABLE  I 


Author(s) 

Marcum 

Marcum 

Raso  and 
Woolf 

George  and 
Lavagnino 

Shelton  and 
Keith 

Raso 

Ce'lnik 

Eamon 


nte  Carlo  Calculations  In  Variable  Density  Air 


Year 

Ref. 

No. 

Source  Height  (km) 

Source 

1963 

2 

30.48,  60.96 

14  MeV  neutron 

1963 

3 

60.96 

51.82 

3 MeV  gamre 
fission  gamma 

1966 

4 

20.00 

1.25  MeV  gamma 

1966 

5,  6 

9.14,  15.24 

12  monoenergetic 
neutron  sources: 
14.0-0.004  MeV 

1969 

n 

i 

25.91,  33.53 

8 monoenergetic 
neutron  sources: 
14.0-0.03  MeV 

1969 

8 

20,  30,  40 

0.5  and  1.0  MeV 
neutron 

1969 

9 

20  to  75 

7 monoenergetic 
neutron  sources: 
14.1-1.5  MeV 

1976 

10 

20  , 30  , 40  , 60, 
80 

fission  and  thermo- 
nuclear neutron 
sources 

3 


The  application  of  Monte  Carlo  methods  to  this  problem  is  a science 
(some  say  art)  in  its  own  right  and,  unfortunately,  expensive.  For 
example,  Eamori  (Ref.  10)  reports  that  CDC-7600  computer  execution  times 
of  30-40  minutes  are  required  to  obtain  acceptable  statistics  (standard 
deviations  less  than  25%)  for  radiation  environments  from  a fission 
source  at  150  receivers  located  at  mass  ranges*  varying  from  10  to  200 
g/cm2  of  air. 

The  high  computational  costs  associated  with  Monte  Carlo  determina- 
tions coupled  with  the  multi  scenario  nature  of  nuclear  S/V  assessments 
led  to  the  development  and  widespread  use  of  quick-running,  user-oriented 
mass  scaling  codes  like  SMAU6  (Ref.  II),  ATR  (Ref.  12),  and  CDR  (Ref.  13). 
All  of  these  codes  use  mass  integral  scaling  (MIS)  of  infinite  homogeneous 
air  transport  data  to  generate  radiation  environments  which  purport  to 
reflect  the  exponential  nature  of  two-dimensional  (2-D)  air.  Their 
required  input  normally  involves  only  a specification  of  the  source  in 
terms  of  yield  and  spectrum,  source  location,  and  receiver  locations. 

They  are  capable  of  calculating  environments  at  hundreds  of  receiver 
locations  in  a few  seconds. 

Past  attempts  to  evaluate  the  accuracy  of  the  MIS  approximation 
have  been  limited  in  number  and  scope.  In  1963,  Marcum  (Refs.  2 and  3) 
compared  mass  integral  scaled  doses  from  a 14  MeV  neutron  and  a 3 MeV 
gamma  source  to  Monte  Carlo  determined  2-D  air  doses  at  source  altitudes 
of  30.48  and  60.96  km.  He  found  that,  at  these  altitudes,  particles 


*Mass  range  is  the  integral  of  the  density  along  a straight  line  between 
the  source  and  a receiver  point.  In  constant  density,  homogeneous  air  it 
is  simply  the  product  of  the  density  and  the  slant  range  between  the  source 
f.nd  receiver. 
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streaming  out  of  the  top  of  the  exponential  atmosphere  (leakage)  resulted 
in  doses  that  were  as  much  a factor  of  two  lower  thun  the  scaled  doses. 

In  1966,  Raso  and  Woolf  (Ref.  4)  reported  that  MIS  seemed  to  work  pretty 
well  for  a 1.25  MeV  gamma  source  at  an  altitude  of  20  km.  Unfortunately, 
their  evaluation  was  restricted  to  mass  ranges  less  than  80  g/cm?. 

In  1969,  several  evaluations  of  MIS  were  reported.  Raso  (Ref.  8) 
concluded  that  the  MIS  errors  for  a 0.5  MeV  and  a 1.0  MeV  neutron  source 
ac  altitudes  of  20,  30,  and  40  km  were  small.  These  results  are  some- 
what misleading  because  his  evaluation  was  restricted  to  small  mass  ranges 
of  less  than  20  g/cm2.  Celnik  (Ref.  9)  evaluated  MIS  for  a number  of 
monoenergetic  neutron  sources  ranging  in  energy  from  1.5  to  14.1  MeV  and 
at  altitudes  ranging  from  20  to  75  km.  His  results  showed  that  the 
leakage  effect  cited  above  resulted  in  order  of  magnitude  errors  for 
source  altitudes  above  40  km  and  ranges  greater  than  80  g/cm2. 

It  should  be  mentioned  that  the  ATR  code  (Ref.  12)  has  a correction 
factor  that,  purports  to  correct  scaled  doses  for  exponential  air  effects. 
The  correction  algorithm  is  based  on  an  extension  of  French's  (Ref.  14) 
'First-Last  Collision  Model"  which  was  developed  to  predict  the  air/ 
ground  interface  perturbation  on  doscis  from  fast  neutrons.  An  examina- 
tion of  the  algorithm  shows  that  the  ATR  correction  factor  is  a function 
only  of  the  source  and  receiver  altitude  and  is  indeiiendent  of  range. 

In  contrast,  the  results  of  the  evaluations  cited  above  show  quite 
clearly  that  the  correction  factor  is  strongly  dependent  on  range.  The 
ATR  correction  also  assumes  that  the  only  perturbation  Is  the  leakage 
effect  mentioned  above.  Again,  the  results  of  past  evaluations  (Refs.  2, 
3,  7)  have  shown  that  differences  in  mass  distribution  between  the  source 
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and  receiver  in  homogeneous  and  variable  density  air  can  cause  an  addi- 
tional effect  that  is  quite  different  from  leakage.  This  effect,  here- 
after referred  to  as  the  mass  distribution  effect,  can  actually  result 
in  the  2-D  variable  density  air  results  being  higher  than  the  scaled 
results  rather  than  lower  as  leakage  would  predict. 

By  far,  the  most  comprehensive  evaluation  to  date  was  reported  by 
Shelton  and  Keith  (Ref.  7)  in  1969.  They  compared  the  scaled  total 
fluence  from  eight  monoenergetic  neutron  sources  ranging  in  energy  from 
0.03  to  14.1  MeV  to  Monte  Carlo  determined  variable  density  air  flusnces. 
Source  altitudes  of  25.91  and  33.53  km  ard  receiver  altitudes  ranging 
from  18  to  53  km  were  considered.  In  many  cases,  their  evaluation 
extended  to  depths  of  120  g/cm*.  Their  results  again  showed  that  signifi- 
cant differences  (errors  >505)  existed  between  the  scaled  and  actual 
variable  density  air  results. 

In  stannary,  a nunfcer  of  attempts  have  been  made  in  the  past  to 
evaluate  the  accuracy  of  the  MIS  approximation  in  variable  density  air. 
These  evaluations  have  been  limited  in  scope  in  that  only  a few  mono- 
energetic  sources  and  altitudes  were  considered.  MIS  has  not  been 
evaluated  for  spectra  such  as  fission  or  thermonuclear  sources,  for 
source  altitudes  less  than  20  km,  or  for  the  very  Important  secondary 
gamma  environments.  Further,  an  acceptable  MIS  correction  factor 
algorithm  has  not  been  developed. 

Purpose 

The  purpose  of  the  research  reported  here  was  to  conduct  a systematic 
evaluation  of  the  MIS  technique  as  applied  to  the  atmospheric  radiation 
transport  problem.  The  evaluation  is  based  on  determining  the  ratio  of 
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the  2-D  4itRa  variable  density  air  dose  to  the  1-D  homogeneous  air  MIS 
dose  as  a function  cf  mass  range,  for  a particular  source,  dose  response 
function,  and  source/receiver  altitude.  Silicon  and  tissue  dose  ratios 
were  determined  for  neutron  and  secondary  gamma  environments  from  both 
an  unclassified  fission  and  unclassified  thermonuclear  source  for  source 
altitudes  between  5 and  20  km,  for  receiver  altitudes  between  1 and  25 
km,  and  out  to  mass  ranges  of  250  g/cm2. 

The  unique  and  distinguishing  feature  of  this  research  is  that  the 
evaluation  was  conducted  with  a special  form  of  diffusion  theory  which 
Is  defined  later.  The  complete  parametric  evaluation  via  diffusion 
theory  was  conducted  at  a fraction  of  the  cost  ("1/30)  that  would  have 
been  required  by  a comparable  Monte  Carlo  evaluation.  Further,  by 
comparing  diffusion  derived  results  to  available  Monte  Carlo  result  . 
it  Is  shown  that  the  diffusion  evaluation  is  an  adequate  model  for 
predicting  the  MIS  evaluation  ratios. 

This  report  contains  seven  sections  and  six  appendixes.  In  Section 
II,  the  mass  scaling  law  is  described,  and  the  application  of  MIS  in  the 
scaling  codes  is  discussed.  The  method  used  to  evaluate  MIS  is  described 
in  detail  in  Section  III.  The  diffusion  equation  in  variable  density 
air  is  developed  in  Section  IV.  Section  V contains  an  examination  of 
the  validity  of  the  diffusion  calculations.  Comparisons  of  diffusion 
derived  radiation  environments  and  MIS  evaluation  ratios  are  discussed 
here.  The  actual  evaluation  uf  mass  integral  scaling  is  presented  in 
Section  VI.  Finally,  a summary  and  concluding  remarks  are  presented 
in  Section  VII. 
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I I . Mass  Integral  Scaling 


The  Mass  Scaling  Law 

The  m,  -S  scaling  law  is  widely  used  and  is  rigorously  applicable  to 
the  transport  of  radiation  in  homogeneous  media.  For  application  to  the 
atmospheric  radiation  transport  problem,  it  can  be  stated  as  follows: 

In  an  infinite  homogeneous  medium  with  an  isotropic  point  source,  the 
4uR2  fluence  (time  integrated  flux)  or  dose  is  a function  only  of  pR,* 
the  mass  per  unit  area  between  the  source  and  a receiver.  In  1956, 

Zerby  (Ref.  15)  presented  a very  simple  derivation  of  the  mass  scaling 
law  using  the  Boltzmann  transport  equation.  As  shown  in  Appendix  A, 
the  scaling  law  can  be  derived  from  the  diffusion  equation  in  a similar 
manner  and  under  the  exact  same  set  of  assumptions. 

To  further  clarify  how  mass  scaling  is  applied  in  homogeneous  air, 
let's  consider  the  following  example.  If  identical  isotropic  point 
sources  are  placed  in  two  infinite  homogeneous  air  media  (medium  A and 
medium  B),  and  if  two  slant  ranges,  and  Rg,  in  the  two  ' edi a are 
relat.ed  by 

pBRB  = PARA 

where  p is  the  density  in  the  respective  media,  then,  as  shown  in 
Appendix  A,  the  energy  dependent  fluences,  F(R,E),  at  ranges  R^  and  Rg 

*In  homogeneous  air,  "pR"  is  the  mass  range  or  areal  density  and  is  given 
by  the  product  of  the  density  and  slant  range  between  the  source  and 
receiver. 


8 


must  be  related  as 

, 4ttR|  F(RB,E)  - 4itR*  f(ra,e)  (2.2) 

It  is  convenient  to  refer  to  the  quantity  on  the  left-  (or  right-)  hand 
side  of  equation  (2.2)  as  the  "47rR2  fluence." 

If  the  atmosphere  were  homogeneous,  the  mass  scaling  law  could  be 
rigorously  applied,  without  error,  to  quickly  compute  radiation  environ- 
ments under  a wide  variety  of  scenarios.  Unfortunately,  this  is  not  the 
case;  the  density  in  the  U.S.  Standard  Atmosphere,  as  shown  in  Figure  1, 
is  not  constant  but  decreases  in  an  exponential-like  fashion  with  altitude. 
In  the  next  subsection,  it  is  shown  how  the  mass  scaling  law  is  extended 
to  variable  density  air  and  becomes  an  approximation  whose  accuracy 
varies  and  depends  on  a number  of  parameters. 


Extension  of  the  Mass  Scaling  Law  to  Variable  Density  Air 

Since  the  density  is  not  constant  in  variable  density  air,  defined 
here  to  be  the  U.S.  Standard  Atmosphere  (Ref.  1),  the  application  of  mass 
scaling  here  requires  a careful  definition  of  mass  range.  The  mass  range 
in  this  atmosphere  is  defined  to  be  the  mass  integral,  <pR>,  which  can  be 
written  as 


<pBrB> 


(2.3) 


and  is  also  sometimes  referred  to  as  the  areal  density. 

To  show  how  the  mass  scaling  law  is  applied  in  variable  density  air, 
assume  that  medium  A is  homogeneous  air  and  that  medium  B is  the  U.S. 
Standard  Atmosphere.  The  mass  integral  scaling  (MIS)  approximation 
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assumes  that  If  the  mass  range  between  a source  and  receiver  in  variable 
density  air  is  equal  to  the  mass  range  in  homogeneous  air,  then  the  4irR2 
fluences  will  be  the  same,  i.e.,  if 

<pBRB>  = pARA 

then 

4ttR2  F(RB,E)  = 4ttR2  f(ra,e)  (2.5) 

So,  the  fundamental  assumption  made  in  MIS  codes  is  that  homogeneous 
air  transport  results  can  be  mass  integral  scaled  to  generate  radiation 
environments  which  purport  to  reflect  the  near-exponential  nature  of  the 
U.S.  Standard  Atmosphere.  All  of  these  codes  have  an  infinite  homogeneous 
air  transport  data  base  like  the  one  presented  by  Burgio  in  reference  16. 
The  data  base  is  stored  as  4irR2  fluence  as  a unction  of  mass  range  (areal 
density).  To  calculate  the  approximate  2-D  atmospheric  dose,  a user 
specifies  the  source  spectrum,  dose  response  function,  source  height, 
receiver  height,  and  slant  range.  The  code  then  evaluates  the  mass 
integral  and  uses  its  data  base  to  compute  the  1-D  MIS  4nR2  dose  at  this 
mass  range.  Finally,  the  dose  is  found  by  simply  dividing  by  4ttR2. 

Deficiencies  of  Mass  Scaling  in  Variable  Density  Air 

Since  the  mass  scaling  law  was  derived  under  the  assumption  of  an 
infinite  homogeneous  media,  it  should  be  clear  that  it  becomes  an  approxi- 
mation when  extended  to  variable  density  air.  The  most  ob”'ous  deficiency 
of  MIS  when  applied  in  this  medium  is  its  inability  to  account  for  the 
streaming  of  particles  out  the  top  of  the  exponential -like  atmosphere 
(leakage).  In  homogeneous  air,  the  medium  is  infinite  in  that  every  point 
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is  theoretically  surrounded  by  an  infinite  amount  of  mass.  Therefore, 
it  is  not  possible  for  particles  to  "leak."  In  the  U.S.  Standard  Atmos- 
phere, the  density  decreases  in  an  exponential-like  fashion  with  altitude. 
Figure  2,  which  depicts  the  integral  with  respect  to  altitude  of  the 
density  function  given  in  Figure  1,  shows  that  the  mass  of  air  above  an 
altitude  z also  decreases  exponentially  with  altitude.  At  altitudes 
above  40  km,  there  are  only  a few  grams/cm2  of  air  above,  and  particles 
directed  upward  at  these  altitudes  can  and  do  stream  out  the  top  of  the 
atmosphere.  Since  mass  integral  scaling  cannot  account  for  leakage,  the 
scaled  doses  must  be  higher  than  the  actual  2-D  variable  density  air 
doses.  The  severity  of  the  leakage  effect  obviously  increases  with  alti- 
tude. 

Another  more  subtle  deficiency  of  MIS  is  due  to  the  differences  in 
mass  distribution  in  homogeneous  air  and  in  the  U.S.  Standard  Atmosphere. 

In  homogeneous  air,  which  is  used  to  construct  the  scaling  data  base,  the 
mass  is  uniformly  distributed  between  the  source  and  all  receiver  points. 

In  contrast,  the  mass  is  nearly  exponentially  distributed  in  the  vertical 
direction  in  the  U.S.  Standard  Atmosphere.  In  this  medium,  a substantial 
number  of  particles  can  reach  receivers  after  traveling  great  distances 
in  the  less  dense  air  above  the  source.  Although  their  paths  are  much 
greater  than  the  simple  slant  range  between  the  source  and  these  receivers, 
the  particles  encounter  much  less  air  in  traversing  these  distances  than 
predicted  by  the  simple  mass  integral.  Therefore,  they  suffer  less 
attenuation.  The  net  effect  of  the  mass  distribution  difference  is  to 
enhance  the  2-D  variable  density  air  dose  at. large  mass  ranges  and  depress 
it  at  smaller  ranges.  The  mass  distribution  effect  is  clearly  present  in 
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Marcum's  data  (Refs.  2,  3)  and  has  been  discussed  previously  by  Shelton 
and  Keith  (Ref.  7). 

Finally,  MIS  is  clearly  inadequate  for  accurately  predicting  doses 

t 

from  low  altitude  bursts  and  for  receivers  near  the  air/ground  interface. 
This,  of  course,  is  due  to  the  change  of  medium  and  drastic  discontinuity 
of  density  that  occurs  at  the  air/ground  interface.  Straker's  (Ref.  17) 
analysis  showed  that  the  ground  acts  as  a reflector  at  ranges  close  to 
the  source  and  as  an  absorber  at  ranges  far  from  the-  source.  Therefore, 
the  scaled  doses  are  low  in-close  and  high  at  greater  ranges.  Recently, 
Huszar  (Ref.  18:38)  compared  some  actual  2-D  air/ground  transport  calcu- 
lations (Ref.  19)  to  scaled  results.  He  reports  that  if  the  source/ 
receiver  altitude  is  greater  than  1000  meters  and  the  corresponding 
receiver/source  altitude  is  greater  than  200  meters,  then  the  air/ground 
perturbation  of  the  scaled  dose  is  negligible. 

In  summary,  the  effects  of  leakage,  mass  distribution,  and  the  air/ 
ground  interface  can  all  produce  perturbations  in  the  mass  integral 
scaled  dose  in  variable  density  air.  As  stated  at  the  outset,  the  purpose 
of  the  research  reported  here  was  to  conduct  a systematic  evaluation  of 
the  MIS  technique  as  applied  to  the  atmospheric  radiation  transport 
problem.  Air/ground  perturbations,  which  have  recently  been  examined  in 
detail  (Refs.  18,  19),  were  not  evaluated  in  the  effort  reported  here. 
Results  will  be  presented  later  which  clearly  define  the  error  due  to 
leakage  and  mass  distribution  over  various  spatial  domains  in  the  1-25  km 
altitude  region.  In  the  next  section,  the  method  of  evaluation  used  in 
this  study  is  defined. 
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III.  Method  of  Evaluation 


The  K Factor 

Many  nuclear  S/V  assessments  do  not  require  such  details  as  the 
angular  or  energy  dependent  fluence  and  are  conducted  using  scalar  quan- 
tities like  dose  which  is  defined  by 

Dose(r)  = ]£  Fs(r)D§  (3.1) 

g=i  R 


where 

F9(r)  * the  angle  integrated  scalar  fluence  at  a distance  of  |r| 
from  the  source  for  energy  group  g of  a multi  group  energy 
set  which  contains  G groups  (particles/cm2). 
d|  * the  experimentally  determined  dose  response  function  for 
group  g in  units  of  dose  (e.g.,  rads  silicon  or  tissue)/ 
unit  scalar  fluence. 

Therefore,  in  this  study  mass  integral  scaling  was  evaluated  by  deter- 
mining K factors  which  were  defined  as 


K factor  * 2-0  4itR2  VARIAplE  DENSITY  AIR  DOSE  2j 

1-D  4ttR2  MIS  DOSE 

Since  the  4ttR2  variable  density  air  dose  is  a function  of  the  source 
height,  receiver  altitude,  mass  range,  source  type,  and  response  function, 
the  K factor  must  also  be  dependent  on  these  parameters.  This  evaluation 
method  is  convenient  because  the  K factor  can  be  used  not  only  as  a 
measure  of  the  error  due  to  MIS  but  also  as  a simple  multiplicative 
correction  factor  for  the  scaled  doses.  The  error  due  to  MIS,  which  is 
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defined  to  be  the  absolute  value  of  the  difference  between  the  2-D  and 
1-D  4nR2  doses  divided  by  the  2-D  4ttR2  dose,  is  related  to  the  K factor  by 


K factor  - 1 
K factor 


(3.3) 


Calculation  Parameters 

In  this  study,  K factors  were  determined  for  neutrons  and  secondary 
gammas,  for  a fission  and  a thermonuclear  source,  for  source  altitudes 
between  5 and  29  km,  for  receiver  altitudes  between  1 and  25  km,  for  mass 
ranges  out  to  250  g/cm2,  and  for  silicon  and  tissue  doses.  The  calcula- 
tions reported  here  were  made  using  the  DLC-31,  37/21  coupled  neutron  and 
gamma,  multi  group  energy  structure  (Ref.  20)  shown  in  Table  II,  the 
silicon  and  tissue  response  functions  shown  in  Tables  III  and  IV,  and  the 
unclassified  fission  and  thermonuclear  sources  shown  in  Table  V.  The 
DLC-31  fission  weighted  multi  group  cross  section  set  described  in  reference 
16  was  used.  Macroscopic  cross  sections  were  obtained  under  the  assump- 
tion that  air  at  a density  of  1.11  mg/cm3  is  composed  of  79%  nitrogen  and 
21%  oxygen  by  volume  (3.6609  x 1019  N atoms/cm3  and  9.7316  x 1018  0 atoms/ 
cm3). 

As  pointed  out  earlier,  the  unique  and  distinguishing  feature  of  this 
research  is  that  the  computation  of  K factors  and  MIS  evaluation  were  con- 
ducted with  a special  form  of  diffusion  theory.  This  approach  was  chosen 
because  an  insufficient  number  of  Monte  Carlo  calculations  existed  (Table  I). 
Further,  the  cost  to  make  the  number  of  Monte  Carlo  calculations  that  might 
be  required  to  adequately  account  for  the  sensitivity  of  the  K factor  to 
source/response  type,  source/receiver  altitude,  and  range  was  thought  to 
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be  prohibitive.  The  special  form  of  diffusion  theory  Is  developed  In  the 
next  section,  and  the  adequacy  of  Its  use  In  the  MIS  evaluation  Is  examined 
In  Section  V. 


TABLE  IIa 

The  DLC-31  Multi  group  Energy  Structure  and  Group  Mean  Free  Paths 

Upper  Energy  Mean  Free  Path0  Upper  Energy  Mean  Free  Path 

Group  (MeV)  (g/cm2)  Group  (MeV)  (q/cm2) 


Neutron  Groups 


1 

1. 9640+1 b 

15.2 

20 

1.8268 

10.9 

2 

1 .6905+1 

15.1 

21 

1 . 1080 

10.7 

3 

1.4918+1 

15.0 

22 

5.5023-1 

6.3 

4 

1.4191+1 

15.2 

23 

1.5764-1 

6.1 

5 

1.3840+1 

15.2 

24 

1.1109-1 

5.2 

6 

1.2040+1 

15.3 

25 

5.2475-2 

4.3 

7 

1.2214+1 

16.1 

26 

2.4788-2 

4.0 

8 

1.  352+1 

17.1 

27 

2.1875-1 

3.7 

9 

1.0000+1 

18.6 

28 

1 .0333-2 

3.5 

10 

9.0484 

19.0 

29 

3.3546-3 

3.2 

11 

8.1873 

17.0 

30 

1.2341-3 

3.1 

12 

7.4082 

18.8 

31 

5.8295-4 

2.9 

13 

6.3763 

17.3 

32 

1.0130-4 

2.8 

14 

4.9659 

21.6 

33 

2.9023-5 

2.8 

15 

4.7237 

13.3 

34 

1 .0677-5 

2.7 

16 

4.0657 

12.5 

35 

3.0590-6 

2.7 

17 

3.0119 

17.5 

36 

1.1254-6 

2.7 

18 

2.3852 

19.9 

37 

4.1400-7 

2.4 

19 

2.3069 

14.8 

aFrom  reference  16. 

bRead  1.9640+1  as  1.9640  x 10*. 

cThe  mean  free  path  In  mass  units  is  Independent  of  altitude  and  computed 
by  multiplying  the  mean  free  path  (1/Zt)  at  altitude  z by  the  density  at 
altitude  z. 


TABLE  II  (continued) 

The  DLC-31  Nultl group  Energy  Structure  and  Group  Mean  Free  Paths 


Group 

Upper  Energy  1 
(MeV) . 

Mean  Free  1 

- Ja/gg2 

Ganna  Groups 

38/1 

1.4000+1 

52.3 

39/2 

1 .0000+1 

47.4 

40/3 

8.0000 

44.2 

41/4 

7.0000 

41.4 

42/5 

6.0000 

38.2 

43/6 

5.0000 

34.5 

44/7 

4.0000 

30.3 

45/8 

3.0000 

26.7 

46/9 

2.5000 

24.0 

47/10 

2.0000 

20.9 

48/11 

1.5000 

17.5 

49/12 

1.0000 

14.5 

50/13 

7.0000-1 

12.1 

51/14 

4.50000-1 

10.2 

52/15 

3.0000-1 

8.4 

53/16 

1.5000-1 

7.1 

54/17 

1.000C  -1 

6.4 

55/18 

7.0000-2 

5.6 

56/19 

4.5000-2 

4.4 

57/20 

3.0000-2 

2.6 

58/21 

2.0000-2 

0.7 
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TABLE  III 


Neutron  Dose  Response  Functions 


Group 

Silicon  Response® 
Rads/(n/cm2) 

Tissue  Response 
Rads/(n/cml) 

1 

1.9106E-09 

8.67241E-09 

2 

1.7792E-09 

7.41 901 E-09 

3 

1.6818E-09 

6.811 54E-09 

4 

1.6231 E-09 

6. 54471 E-09 

5 

1.5144E-09 

6.14729E-09 

6 

1 .3851 E -09 

5.95478E-09 

7 

1.2370E-09 

5.89361E-09 

8 

1.0530E-09 

5.55077E-09 

9 

8.7897E-10 

5.28824E-09 

10 

7.9629E-10 

5. 04731 E-09 

11 

7.8141E-10 

5.00452E-09 

12 

4.7092E-10 

4.75953E- 09 

13 

2.1394E-10 

4.48311E-09 

14 

1.8267E-10 

4. 25311 E-09 

15 

1 .41 95E-1 0 

4.1711 0E-09 

16 

1.0582E-10 

3.97837E-09 

17 

1.0006E-10 

3.39051 E-09 

18 

8.2995E-11 

3.13766E-09 

19 

9.4778E-U 

3.03449E-09 

20 

6.5328E-U 

2.63932E-09 

21 

4.9785E-11 

2.05699F-09 

22 

3.1515E-11 

1.33305E-09 

23 

1 . 7897E-1 2 

7.62280E-10 

24 

2.8022E-12 

5.48897E-10 

25 

1.2327E-12 

3.11 644E-1 0 

26 

7.9084E-13 

2.07393E-10 

27 

5.8930E-13 

1.46624E-10 

28 

2.9804E-13 

6.61426E-11 

29 

1.0498E-13 

2.27582E-11 

30 

4.3305E-14 

9.13150E-12 

31 

1.4421E-14 

3.66325E-12 

32 

4.5895E-15 

1.17590E-12 

33 

3.9377E-15 

1.1 0948E-1 2 

34 

5.6286E-15 

1 .61171 E- 1 2 

35 

9.4023E-15 

2 . 741 57E-1 2 

36 

1.5380E-14 

4.45699E-12 

37 

7.4244E-13 

1.1 2382E-1 1 

*Prov1ded  by  Dave  Bartlne  of  ORNL. 

‘Based  on  rebanding  of  data  In  reference  21  by  H.  M.  Murphy. 
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T4GLE  IV 


6mm  Dose  Response  Functions 


Silicon  Response5 

Tissue  Response 

Group 

Rad/ (6/ cm2) 

Rad/ (6/ cm*) 

1 

3.41838E-09 

2.7431 1 E-09 

2 

2.57121E-09 

2.25643E-09 

3 

2.16121E-09 

1.98395E-09 

4 

1 -89915E-09 

1.79216E-09 

5 

] .63669E-09 

1.59277E-09 

6 

1.38348E-09 

K38972E-09 

7 

1 . 13348E-09 

1.18027E-09 

8 

9.43338E-10 

1.00980E-09 

9 

8.20339E-10 

8.83204E-10 

10 

6.83433E-10 

7. 4281 0E- 10 

11 

5.?8457E-10 

5.80300E-10 

12 

3.85046E-10 

4.23930E-10 

13 

2.71224E-10 

2.96950E-10 

14 

1.77718E-10 

1.92834E-10 

15 

1.04590E-10 

1 .07596E-10 

16 

6.85104E-1 1 

4.93830E-11 

17 

7.98544E-U 

3.43145E-1 1 

18 

1.45429C-10 

2.94789E-11 

19 

3.44112E-10 

4.37499E-11 

20 

8.26793E-10 

9.66468E-11 

21 

2.64927E-09 

3.25042E-10 

aBased  on  rebanding  of  data  given  in  reference  22  by  H.  M.  Murphy. 
bBased  on  rebanding  of  data  given  in  reference  23  by  H.  M.  Murphy- 
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TABLE  V 


Neutron  Sources3 
Group  Thermonuclear 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

Total 

®From  reference  24 


0.0 

0.0 

1.887 i 4-02* 

9.34254-03 

2.66169-02 

1.66622-02 

1.68678-02 

1.23974-02 

7.48258-03 

6.82320-03 

6.77521-03 

1.03201-02 

1.80706-02 

3.61700-03 

1.24302-02 

2.60380-02 

2.37305-02 

3.74662-03 

2.56418-02 

6.44472-02 

8.84954-02 

9.13765-02 

1.16335-02 

1.10777-01 

5.40049-02 

5.68196-03 

9.26377-02 

1.16267-01 

7.38166-02 

2.32454-02 

2,02810-02 

1.90145-03 

0.0 

0,0 

0.0 

0,0 

0.0 

1.00000 

♦Read  1.88714-02  as 


Fission 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3.83946-03 

3.50150-03 

5.38923-03 
7.34895-03 
1 .83674-02 
3.24954-03 
8.46880-03 

5.50022-02 

3.24354-02 

1.05785-02 

9.72362-02 

1.46769-01 

2.15670-01 

1.50179-01 

1.92973-02 

1.20983-01 

5.72919-02 

5.99914-03 

2.39968-02 

1.43980-02 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.00000 
.88714  x 10"2 
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IV.  Atmospheric  Diffusion  Theory 


This  section  contains  a development  of  the  atmospheric  diffusion 
theory  used  in  this  research.  There  are  twc  features  of  the  theory  which 
are  unique  and  deserve  special  mention.  First,  the  1962  U.S.  Standard 
Atmosphere  (Ref.  1)  was  incorporated  directly  into  the  diffusion  equation. 
Secondly,  a nunorthogonal  coordinate  system,  which  expands  exponentially 
with  altitude,  was  developed  and  used.  As  will  be  shown  below,  this 
coordinate  system  is  essential  to  the  numerical  solution  of  the  diffusion 
or  transport  equation  in  variable  density  air. 

Two  additional  features,  which  were  included  to  minimize  the  inher- 
ent limitations  of  diffusion  theory,  are  also  worthy  of  mention.  First, 
although  the  group-to-group  scatter  was  assumed  to  be  isotropic  as 
required  by  diffusion  theory,  the  in-group  scatter  was  treated  as  aniso- 
tropic (limited  to  linear  angular  variations)  by  applying  a transport 
correction  to  the  diffusion  coefficient.  Secondly,  diffusion  theory 
restricts  the  angular  dependence  of  the  fluenci  to  linear  variation  with 
the  cosine  of  the  angle,  0 , which  is  defined  by  arcosine  (r  • v),  where 
r is  the  unit  position  vector  and  v is  the  unit  paiticle  velocity  vector 
(Ref.  25:234).  Since  the  uncollided  particles  (direct  fluence)  emanating 
from  an  isotropic  point  source  are  all  directed  radially  outward,  their 
angular  variation  must  be  represented  by  a dirac  delta  function  in  cos  eo 
which,  of  course,  is  poorly  treated  by  the  diffusion  theory  assumption 
of  linearity  with  cos  0Q.  This  limitation  was  circumvented  by  using  an 
analytical  first  collided  source.  With  this  source,  diffusion  theory  is 
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applied  only  to  particles  which  have  scattered  at  least  one  time*  and  the 
direct  or  uncoil ided  portion  of  the  fluence  is  computed  analytically. 

The  Multi  group  Diffusion  Equation 

The  multi  group  diffusion  equation  in  general  coordinates,  xx,  x2, 
x3,  can  be  written  as 

7 • (d9(  Xi,x2,x3)  VFg(xi,x2,x3  >)  - er<  Xi,x2,x3)  Fs(xx,x2,x3) 

■”  -S9( Xx,x2,x3) 


- X)  s?  >9(xiux2,x3)  F9  (xi,x2,x3) 
g*=l  s 

(4.1) 

where 

g = the  energy  group  number  which  varies  from  1 for  the  highest 
energy  to  G for  the  lowest  energy 
D9  » the  group  g diffusion  coefficient  in  units  of  cm 
F9  = the  angle  integrated  group  g scalar  fluence  in  units  of 
particles/cm2 

= the  macroscopic  removal  cross  section  for  group  g in  units 


of  cm”1 


Z9  = the  macroscopic  total  interaction  cross  section  for  group 
g in  units  of  cm"1 

I9  ^ = the  macroscopic  cross  section  for  scatter  from  group  g1  to 
group  g in  units  of  cm"1 
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S9  = the  source  of  particles  for  group  g in  units  of  particles/ 
cm* 

It  should  be  noted  that  the  assumption  of  no  up-scatter  in  energy  has 
been  made  in  writing  equation  (4.1). 


A Nonorthogonal  Coordinate  System 

McLaren  (Ref.  26)  was  the  first  to  recognize  that  a nonorthogonal 
coordinate  system  was  essential  to  the  solution  of  the  diffusion  or 


transport  equation  in  variable  density  air.  This  requirement  arises 
from  the  exponential-like  density  variation  of  the  U.S.  Standard  Atmos- 
phere (Figure  1).  Since  the  density  decreases  exponentially  with  alti- 
tude, the  particle  mean  free  paths*  increase  exponentially  with  altitude. 
As  shown  below,  this  phenomenon  imposes  a unique  problem  on  the  numerical 
solution  of  equation  (4.1). 

Consider  Figure  3a  which  shows  a generalized  set  of  mesh  lines  in  a 
right  circular  cylindrical  coordinate  system  which  might  be  used  in  the 
numerical  solution  of  equation  (4.1)  in  2-D  variable  density  air.  The 


vertical  mesh  spacing  is  set  at  some  arbitrary  fraction  of  a mean  free 
path  (e.g.,  three  lines  per  mean  free  path)  and  increases  with  altitude. 


because  the  mean  free  path  increases  with  altitude.  Although  the  hori- 
zontal meshing  is  shown  to  be  constant  in  Figure  3a,  classical  finite 
differencing  methods  do  permit  the  radial  mesh  to  vary  in  the  radial 


♦The  mean  free  path  can  be  conveniently  expressed  in  altitude  independent 
mass  units  by  multiplying  the  mean  free  path  at  altitude  z by  the  density 
at  altitude  z.  As  shown  earlier  in  Table  II,  the  DLC-31  neutron  groups 
have  mean  free  paths  which  vary  from  2 to  20  g/cm2  while  a 1-50  g/cm2 
variation  exists  for  the  gamma  groups. 


direction.  However,  variation  of  the  radial  mesh  with  altitude,  which 
Is  clearly  required  if  a constant  mean  free  path  separation  between  mesh 
points -is  to  be  maintained,  is  not  permitted.  In  practice,  the  radial 
mesh  spacing  would  be  fixed  by  the  bottom  row  where  the  air  is  most  dense. 
This  would  result  in  a radial  mesh  which  would  become  finer  in  a mean 
free  path  sense  with  increasing  altitude.  This  continuously  decreasing  mesh 
(in  mean  free  paths)  introduces  numerical  instabilities  into  the  solution. 

This  problem  can  be  overcome  by  using  a mesh  which  expands  with 
altitude  as  shown  in  Figure  3b.  The  penalty  of  this  approach  is  non- 
orthogonal  coordinates.  (Note  that  the  mesh  lines  in  Figure  3b  no 
longer  intersect  at  right  angles.)  With  such  a coordinate  system,  it 
is  possible  to  maintain  a constant  mean  free  path  radial  mesh  spacing 
by  simply  allowing  the  coordinate  system  to  expand  at  a rate  which 
approximates  the  exponential-like  density  variation  of  the  U.S.  Standard 
Atmosphere.  Figure  4,  which  shows  the  coaltitude  range  to  100  g/cm2  as 
a function  of  altitude  in  the  U.S.  Standard  Atmosphere,  indicates  the 
rate  of  expansion  required. 

Finally,  it  is  important  to  note  that  an  additional  advantage  is 
realized  by  the  introduction  of  the  nonorthogonal  coordinates.  Specif- 
ically, solutions  can  be  obtained  over  a much  larger  area  without 
increasing  the  number  of  mesh  points.  This  is  significant  because  this 
Increase  in  solution  domain  does  not  require  any  additional  computer 
memory  or  computation  time. 

The  coordinate  system  developed  for  this  study  is  similar  to  the 
one  proposed  by  McLaren.  It  is  shown  in  Figure  5 and  defined  by  the 
following  transformation  equations: 


Xi  B 

r 

(4.2) 

x2  = 

<t> 

(4.3) 

Xj  = 

z 

(4.4) 

X * 

r ez/,H  cos  4> 

(4.5) 

y * 

r eZ/^H  sin  <j> 

(4.6) 

z = 

z 

(4.7) 

where  x,  y,  z are  standard  rectangular  coordinates  and  M is  the  mesh 
expansion  parameter. 

The  Diffusion  Equation  in  the  U.S.  Standard  Atmosphere 

A simple  and  direct  approach  for  incorporating  the  'J.S.  Standard 
Atmosphere  into  the  diffusion  equation  is  to  express  all  cross  sections, 
which  are  functions  of  altitude  only  in  this  medium,  in  terms  of  sea 
level  values.  If  f(z)  is  the  ratio  of  the  density  at  altitude  z to  the 
density  at  sea  level  (i.e.t  f(z)  = p(z)/pQ),  then 

l|(z)  = f(z)  (4.8) 

J*'**{z)  = f(z)  (4.9) 

Dg(z)  = Dg/f(z)  (4.10) 

whet  Eg  and  Dg  are  all  defined  at  sea  level  (pQ  = 1.225  x 1C“3 
g/cms).  As  pointed  out  at  the  beginning  of  this  section,  the  1962  U.S. 
Standard  Atmosphere  (Ref.  1)  was  adopted  as  the  atmospheric  model  in 
this  stu4y.  An  algorithm,  which  uses  reference  l's  equations  (ideal  gas 
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law  and  hydrostatic  equation)  and  defining  atmospheric  parameter  base 
values,  was  incorporated  into  the  AIRDIF  code  to  generate  the  density 
ratio,  f(z),  at  any  altitude  z. 

An  examination  of  equation  (4.10)  and  the  spatial  leakage  term  in 
equation  (4.1)  shows  that  an  additional  parameter  is  required  in  defining 
the  atmosphere.  Specifically,  the  gradient  of  pQ/p(z)  ^g(z)  = ^-(l/f(z)^ 
as  a function  of  altitude  is  also  required.  In  this  study,  a model  based 
on  the  assumption  of  local  exponential  density  variation  in  some  small 
altitude  increment  A z about  z was  developed  to  generate  the  value  of  g(z) 
at  any  altitude  z.  Under  this  assumption,  it  can  be  easily  shown  that 

g(z)  - (t  /p(z))  = (4.11) 

dz  ' 0 ' 2iz  f ( z) 

where 

Zj  = z + Az  (4..12) 

z2  = z - Az  (4.13) 

Figure  6 shows  how  g(z)  varies  with  altitude. 

Since  the  density  and  hence  all  cross  sections  are  functions  of 
altitude  only,  symmetry  with  respect  to  the  angle,  $,  exists  and  the  multi - 
group  diffusion  equation  given  in  equation  (4.1)  can  be  written  in  the 
U.S.  Standard  Atmosphere  and  in  the  nonorthogonal  coordinate  system  as 

V • Qi-VFg(r,z)j  - Zg  f(z)F9(r,z) 

= -S9(r,z)  - ^ I9'"9  f(z)Fg‘(r,z)  (4.14) 

g'*i 
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tfK*-  - 


The  generalized  form  of  the  spatial  leakage  term,  V • (d(z)  VF(r,z)), 
Is  given  by 

7 . (d(z)  VF(r,z))=  ^D(xj)  /g  gKP  ^ (4.15) 

KP 

where  g are  the  elements  of  the  contravariant  metric  tensor  associated 
with  the  coordinate  transformation  defined  by  equations  (4.2)  through 
(4.7),  and  g is  the  determinant  of  the  metric  tensor,  g^p,  associated 
with  that  transformation. 

Under  the  coordinate  transformation  defined  by  equations  (4.2) 
through  (4.7),  the  elements  of  the  metric  tensor  are 


9..  ■ 

(4.16) 

9l2  ■ 921  s 0 

(4.17) 

_ _ _ _ r 2z/H 

9j  j - 93i  " [f  e 

(4.18) 

n _ „2  2z/H 
g22  -re 

(4.19) 

g2s  * gs2  ® o 

(4.20) 

r2  p2z/H 

9)1  ■ ^ + 1 

(4.21) 

H2 

9 - r2  e4z/H 

(4.22) 

KP 

The  elements  of  the  contravariant  metric  tensor,  g , are 


g11  * 

/ri+e-2z/H\ 

(4.23) 

\H2  / 

g12  ■ 

g21  - 0 

(4.24) 
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(4.25) 

(4.26) 

(4.27) 
(4.23) 


When  equation  (4.15)  is  expanded  using  the  expressions  just  given,  it  can 
be  shown  that  the  multi  group  spatial  leakage  term,  V • (d9(z)  VFq),  can 
be  reduced  to 


(4.29) 

where  superscript  g now  again  is  the  energy  group  designation  and  g(z) 
is  the  gradient  defined  by  equation  (4.11). 

To  permit  linear  isotropy  in  the  ia-group  scatter,  a transport 
correction  was  used  in  defining  the  diffusion  coefficient.  Specifically, 


D9  = 1/3  E9r 


(4.30) 


where 


r9 

H-R 


SjR 


the  transport  cross  section  at  sea  level  for  group  a in 

units  of  cm"1 

E?  - e?-9 
! S 1 
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«*  Pi  coefficient  In  the  Legendre  expansion  of  the  angular 
dependence  of  the  scattering  cross  section  for  group  g 
In-group  scatter 

For  the  reasons  given  at  the  beginning  of  this  section  and  also  to 
avoid  the  numerical  problems  associated  with  using  a point  source,  an 
analytical  first  collided  source  was  used  for  S?(r»z)  in  equation  (4.1  4). 
Speci fi cally, 

S9(r,z)  = £ ^ f(z)  Fj  (r,z)  (4.31) 

g1*!  s v 


where 

Fy  (rfz)  * the  uncollided  fluence  in  group  g' 


(4.32) 


and  where 

Si  * the  number  of  group  g*  particles  emanating  from  an  isotropic 

point  source  at  r * 0 and  z * zs 

R * the  slant  range  between  the  source  at  altitude  zs  and 

receiver  at  altitude  z in  units  of  cm  * ((zs  - z)2  + r2 
62z/h) */2 

Ej  * the  total  macroscopic  cross  section  at  sea  level  (z  * 0) 

<pR>  « mass  range  or  areal  density  defined  by  equation  (2.3) 

It  is  important  to  note  that  the  use  of  the  first  collided  source 
in  equation  (4.14)  means  that  only  the  collided  fluence  is  actually 
found  by  diffusion  theory.  The  total  fluence  is  computed  by  simply 
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adding  the  uncollided  and  collided  fluences,  i.e.. 


, Fj(r,z)  = F9(r,z)  + F9(r,z)  (4.33) 

where 

Fj(r,z)  = total  group  g scalar  fluence 

F9(r,z)  * collided  fluence  found  by  numerically  solving  equation 
(4.14) 

F9(r,z)  ® uncollided  fluence  found  by  analytically  evaluating 
equation  (4.32) 

The  specification  of  boundary  conditions  completes  the  defintion  of 
atmospheric  diffusion  theory.  On  the  sides  of  the  expanding  cylinder 
(r  ■ max)  zero  fluence  boundary  conditions  were  applied.  This  condi- 
tion, of  course,  is  artificial  but  does  not  introduce  significant  error 
as  long  as  the  boundary  is  placed  several  mean  free  paths  (20-30  g/cm2) 
from  the  domain  of  Interest.  At  the  top  and  bottom  of  the  expanding 
cylinder,  zero  return  (partial)  current  conditions  were  applied;  for 
example,  at  z e zmax,  the  condition 

J9(r,z)|  * 0 (4.34) 

|zmax 


was  applied.  In  the  expanding  coordinate  system  and  U.S.  Standard 
Atmosphere,  this  condition  becomes 


Fg(r,z)  , D9  /r2c2z/H  + X'*  / r + 9F^\ 
2f(z)  \H2  / \ H5r  3z  / 


= 0 (4.35) 
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z=zmax 


The  bottom  condition  is  again  artificial  and  is  equivalent  to  assuming 
that  the  air  below  zmin  acts  as  a perfect  absorber.  The  top  condition, 
however,  is  a good  approximation  of  the  true  vacuum  condition  (which 
can't  be  applied  in  diffusion  theory)  as  long  as  zmax  is  set  at  an 
altitude  where  there  is  not  much  air  above  (30  km  in  this  study). 

An  important  characteristic  of  the  equations  just  derived  should 
be  mentioned.  Under  the  conditions, 


H ->■  “ 

(4.36) 

g(z)  = o 

(4.37) 

f(z)  = p(zs)/p0 

(4.38) 

all  of  the  equations  developed  here  reduce  directly  to  a description  of 
the  diffusion  equation  in  normal  cylindrical  coordinates  and  homogeneous, 
constant  density  air.  It  is  this  property  that  enables  the  AIRDIF  code 
described  in  Appendix  B to  compute  diffusion  environments  in  both 
variable  density  and  homogeneous  air. 


Numerical  Solution  of  the  Diffusion  Equation 

Standard  operator  differencing  techniques  (Ref.  27:184)  were  us«d 
to  derive  the  multigroup  finite  differenced  form  of  equation  (4.14).  The 
resulting  nine  point  (five  point  in  homogeneous  air)  difference  equation 
for  interior  points  is 


b-  -i  . . F?  i . , + b.  . « F?  . ■•  + b*.-,  F?,i  . i 

1-1, o-l  1-1»0-1  1,0-1  1,0-1  1+1, 0-1  1+1, 0-1 

+ b-T  . F?  , . + b.  . F?  . + b. , , • fV.-j  . 

"1—1,0  i-1,0  1,0  i ,0  i+1,0  i+1 ,0 


+ bi-l  ,j+l  Fi-1,0+1  + bi.0+l  Fi , j+1  + bi+l,j+l  Fi+l,j+l 


■ -i  djAbaw  XL  fert  rt  .n-, 


A 
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bi-l» j+1 


riAzj-l 


HArAzo(Azj  + Azj-i! 


(4.46) 


bi,j+l 


r (2  + tyjAz,  ,) 

j(Azj  + Azj-l)  V 331 


Azj-1 

Hizj(Azj + Azj-i) 


(4.47) 


bi+l  ,j+l  " ~bi-l,j+l 


(4.48) 


and  where 


q z9'-9s9'  / E9’ 

si  j = fj  — — ~~  exp(-  — <pR>-  j 

9 J Kj  \ Po  ^ 


(4.49) 


and  S9  = the  number  of  particles  in  energy  group  g'  emanating  from  the 
point  source 

The  difference  equations  for  boundary  points  can  be  written  in  a similar 
fashion. 

The  difference  equations  were  written  into  the  matrix  form 


A P - P 


(4.50) 


where 


A is  a block,  tridiagonal  matrix  and  a function  of  the  mesh, 
atmospheric  model,  and  energy  group 
$9  is  the  source  vector  for  group  g and  is  composed  of  the  scatter 
of  both  the  collided  and  uncollided  particles  into  n from 
energy  groups  above  g 


is  the  unknown  scalar  fluence  for  group  g 

The  matrix  equations  were  incorporated  into  the  AIRDIF  code  and 
solved  group-by-group  via  a block-iterative  method  known  as  successive 
line  overrelaxation  (SLOR)  (Ref.  28:199).  As  pointed  out  earlier,  only 
the  collided  fluence  was  found  by  numerically  solving  equation  (4.50). 

The  total  fluence  in  each  group  was  found  by  adding  the  uncollided 
fluence  (equation  (4.32))  and  the  numerically  determined  collided 
fluence. 

The  AIRDIF  Code 

The  AIRDIF  code  is  based  on  the  atmospheric  diffusion  theory  just 
presented.  It  is  capable  of  computing  the  radiation  dose  environments 
in  both  homogeneous  air  and  in  the  U.S.  Standard  Atmosphere  from  atmos- 
pheric nuclear  detonations.  In  addition  to  environment  definition,  it 
is  also  capable  of  directly  computing  the  MIS  K factors.  Appendix  B 
contains  a summary  of  the  detailed  description  of  the  AIRDIF  code  pre- 
sented in  reference  29. 

To  compute  K factors,  an  AIRDIF  homogeneous  air  run  for  the  source 
and  dose  response  of  interest  is  made  first.  The  homogeneous  air  doses 
are  then  least  squares  fit  to  a radiation  transmission  equation  developed 
by  Murphy  (Ref.  16:36), 

T(x)  u exp^A  + Bx  + Cx2  + Dx3/2  + Ex1/2  + Fx1'3  + (J£n(x)) 

(4.51) 

where  T(x)  = 4ttR2  1-D  MIS  dose  at  a mass  range  of  x g/cm2. 

Figures  7 and  8 show  the  fit  and  actual  AIRDIF  data  for  neutron  and 
secondary  gamma  silicon  doses  from  a fission  source  in  homogeneous  air. 
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Figure  8.  Least-Squares  Fit  of  AIRDIF  Secondary  Gamma  Silicon  Doses 
from  a Fission  Source  in  Homogeneous  Air 


These  results  are  typical  and  demonstrate  that  good  fits  to  equation 
(4.51)  can  be  obtained.  Note  that  in  these  and  all  other  figures  to 
follow,  the  symbol  used  for  mass  range  (equation  (2.3))  is  "RHOR." 

The  fit  coefficients  for  both  neutrons  and  gammas  are  included  as 
part  of  the  input  for  the  2-D  variable  density  air  run.  Using  these 
coefficients  and  equation  (4.51),  AIRDIF  computes  the  4ttR2  1-D  MIS  doses 
which  are  required  in  computing  the  K factors. 

It  should  be  emphasized  that  this  parameterization  of  the  homoge- 
neous data  is  essential  to  the  direct  computation  of  K factors  within 
AIRDIF.  It  should  also  be  noted  that  scaling  codes  like  SMAUG  have 
small  central  memory  requirements  because  they  also  utilize  these  kinds 


V.  Validity  of  Diffusion  Calculations 


The  validity  of  the  AIRDIF  diffusion  calculations  was  first  tested 
by  compar'rg  Lae  diffusion  environments  to  discrete  ordinates  transport 
results  in  homogeneous  air  and  to  Monte  Carlo  results  in  variable  density 
air.  Comparisons  were  made  for  neutron  and  secondary  gamma  environments, 
for  silicon  and  tissue  dose  response,  for  a fission  and  thermonuclear 
source,  and  out  to  mass  ranges  of  200  g/cm2.  Next,  to  test  the  adequacy 
of  the  diffusion  model  for  predicting  the  MIS  evaluation  ratios,  compar- 
isons of  K factors  derived  from  diffusion  theory  and  Monte  Carlo  tech- 
niques were  made.  Typical  results  of  these  comparisons  are  presented 
and  discussed  in  this  section. 

Comparisons  of  Radiation  Doses  in  Homogeneous  Air 

In  homogeneous  air,  the  AIRDIF  neutron  and  secondary  gamma  doses 
were  compared  to  Burgio's  (Ref.  16)  discrete  ordinates  results  which 
were  produced  with  the  ANISN  code  (Ref.  30).  Results  were  compared  for 
the  fission  and  thermonuclear  source  shown  in  Table  V and  using  the 
silicon  and  tissue  response  functions  presented  in  Tables  III  and  IV. 

It  is  worth  mentioning  that  the  ANISN  data  are  actually  Murphy's  (Ref. 
16:36)  fits  of  Burgio's  data  to  equation  (4.51).  These  fits  are  in  good 
agreement  with  the  ANISN  calculations  with  maximum  differences  of  about 
5%.  The  AIRDIF  results  were  compared  to  the  ANT3N  fits  rather  than  the 
actual  data  because,  in  the  future,  these  fits  will  become  the  new  data 
base  for  SMAUG.  This  approach  is  consistent  with  the  real  objective  of 
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this  research,  i.e.,  evaluate  the  results  from  mass  scaling  codes. 

Figures  C-l  and  C-2  in  Appendix  C are  typical  of  the  results 
achieved.  For  both  sources  and  response  functions,  the  AIRDIF  diffusion 
results  for  neutrons  are  approximately  10-20%  low  for  mass  ranges  less 
that  20  g/cm2.  Over  the  mass  range  of  20  to  about  140  g/cm2,  the  AIRDIF 
results  are  1-15%  higher.  For  ranges  between  140  and  200  g/cm2,  the 
AIRDIF  results  become  increasingly  lower  than  ANISN  with  maximum  differ- 
ences of  about  30%  occurring  at  200  g/cm2.  The  agreement  in  secondary 
gamma  doses  was  similar.  However,  the  differences  were  somewhat  larger 
with  maximum  errors  of  about  40%  occurring  at  200  g/cm2.  In  general,  the 
agreement  was  better  for  tissue  response  which  emphasizes  the  low  energy 
particle  contribution  (where  the  scatter  is  more  isotropic)  more  than 
the  silicon  response.  Also,  the  agreement  was  a little  better,  as  one 
would  expect,  for  the  fission  source  than  for  the  thermonuclear  source. 

Comparisons  of  Radiation  Doses  in  the  U.S.  Standard  Atmosphere 

Similar  comparisons  were  made  in  the  U.S.  Standard  Atmosphere.  In 
this  medium,  the  comparisons  were  made  to  Eamon's  recent  Monte  Carlo 
calculations  (Ref.  10)  performed  with  the  M0RSAIR  (Ref.  31)  code.  Table 
VI  lists  the  comparisons  which  were  made.  Note  that  AIRDIF  is  about  30 
times  faster  than  M0RSAIR  while  making  determinations  at  10  times  the 
number  of  mesh  or  receiver  points.  This  result  was  not  unexpected  and,  in 
fact,  was  the  primary  reason  why  the  diffusion  approach  was  pursued.  Repre- 
sentative samples  of  Table  Vi's  comparisons  are  presented  in  Appendix  C 
(Figures  C-3  to  C-l 4). 

In  general,  the  M0RSAIR  data  considered  in  these  comparisons  had 
standard  deviations  less  than  10%  for  ranges  less  than  100  g/cm2,  10-20% 
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TABLE  VI 


r 


Summary  of  U.S.  Standard  Atmosphere  Comparisons: 
AIRDIF  versus  MORSAIR 


No. 

Source 

Source 
A1  titude 

Run  Time3 
MORSAIR/ AIRDIF 

Receiver 

Altitude 

Response 

1 

Fission 

5.0 

35. 

5.07 

Silicon 

2 

Fission 

10.0 

30. 

8.84 

Silicon 

3 

Fission 

10.0 

lO.OO 

Silicon 

4 

Fission 

10.0 

12.13 

Silicon 

5 

Fission 

15.0 

32. 

11.30 

Silicon 

6 

Fission 

15.0 

15.29 

Silicon 

7 

Fission 

15.0 

21.66 

Silicon 

8 

Fission 

20.0 

37. 

14.85 

Silicon 

9 

Fission 

20.0 

14.85 

Tissue 

10 

Fission 

20.0 

20.61 

Silicon 

11 

Fission 

20.0 

20.61 

Tissue 

12 

Fission 

20.0 

22.08 

Silicon 

13 

Fission 

20.0 

22.08 

Tissue 

14 

Thermonuclear 

20.0 

29. 

14.85 

Silicon 

15 

Thermonuclear 

20.0 

20.61 

Silicon 

16 

Thermonuclear 

20.0 

22.08 

Silicon 

aRatio  of  MORSAIR  to  AIRDIF  7600  execution  times.  MORSAIR  times  were  for 
50,000  histories  and  150  receivers  (30-40  min).  AIRDIF  times  are  for  on 
the  average  around  1500  mesh  points. 
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for  ranges  between  100  and  150  g/cm2,  and  20-30J  for  ranges  from  150  to 
200  g/cm2.  In  many  cases,  the  AIRDIF/MORSAIR  differences  were  within 
the  statistical  scatter  of  the  Monte  Carlo  data. 

It  is  significant  that  the  agreement  between  diffusion  theory  and 
transport  theory  (Monte  Carlo  in  variable  density  air)  is  about  the 
same  in  the  U.S.  Standard  Atmosphere  as  it  is  in  homogeneous  air.  This 
result  supports  the  validity  of  the  use  of  diffusion  theory  to  compute 
K factors  and  is  clearly  shown  in  the  K factor  comparisons  presented 
next. 


K Factor  Comparisons 

Figures  D-l  through  D-32  in  Appendix  D present  comparisons  of  the 
diffusion  and  transport  (Monte  Carlo)  derived  K factors  for  all  of  the 
detemri nations  listed  in  Table  VI.  For  reasons  wi'ich  are  discussed 


below,  the  agreement  is  better  in  some  cases  than  others.  However,  in 
all  cases  the  K factor  trends,  which  are  indicative  of  the  previously 
discussed  leakage  and  mass  distribution  effects,  are  clearly  the  same. 
The  Monte  Carlo  data  used  in  the  comparisons  were  produced  with  the 


M0RSAIR  code  by  Eamon  of  Kaman  Sciences  Corporation  under  contract 
with  the  Air  Force  Weapons  Laboratory  (Ref.  10).  In  the  odd-numbered 
figures,  D-l  through  D-20,  the  M0RSAIR  and  AIRDIF  neutron  silicon  K 
factors  for  a fission  source  at  source  altitudes  of  5,  10,  15,  and  20  km 


are  compared.  Comparisons  of  the  secondary  gamma  silicon  K factors  are 
made  in  the  even-numbered  figures  in  this  sequence.  Comparisons  are 
also  made  for  tissue  K factors  for  this  source  at  20  km  in  Figures  D-21 
through  D-26.  Finally,  the  silicon  K factors  for  a thermonuclear 


iJ&Ax'l, 


source  at  20  km  are  compared  in  Figures  D-27  through  D-32.  The  sources 
and  response  functions  listed  in  Tables  I I I-V  were  used  in  these  calcu- 
lations. 

The  AIRDIF  and  MORSAIR  neutron  K factors  are  in  excellent  agreement 
with  each  other  for  almost  all  of  the  cases  shown  in  Figures  D-l  through 
D-32.  In  fact,  most  of  the  time  the  AIRDIF  results  are  within  the 
statistics  of  the  MORSAIR  data.  It  is  Important  to  remember  that  the 
AIRDIF  results  were  produced  at  a fraction  of  the  cost  ( ~ 1/30)  required 
by  the  Monte  Carlo  calculations. 

The  agreement  in  the  secondary  gamma  K factors  is  good  but  not 
quite  as  good  as  the  agreement  in  neutron  K factors.  In  general,  the 
AIRDIF  secondary  gamma  K factors  are  lower  than  the  MORSAIR  K factors 
for  ranges  less  than  40  g/cm2  (e.g..  Figure  D-12)  and  higher  for  ranges 
greater  than  140  g/cm2  (e.g..  Figure  D-18).  These  differences  are  due 
to  a combination  of  the  factors  discussed  below. 


Discussion  of  K Factor  Comparisons 

Ideally,  the  K factor  comparisons  should  indicate  only  whether 


where 


K factoryAIRDiF  » K factor) 


MORSAIR 


„ . 2 d)airdif 

K factor/AIRDIF  ' " n\ 

~U/AIRDIF 

* *HmORSAIR  ■ 


(5.1) 


(5.2) 


(5.3) 


(MORSAIR 

and  where  H2-D"  and  "l-D”  are  the  4ttR*  2-D  variable  density  air  and  1-D 
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HIS  doses,,  respectively.  Unfortunately,  the  automation  of  the  K factor 
calculations  In  both  the  AIRDIF  and  MORSAIR  codes  Introduced  several 
additional  factors  which  must  be  considered  In  analyzing  the  results  of 
the  comparisons.  Specifically,  the  quantities  that  are  actually  compared 
In  Figures  D-l  through  D-32  are 

K '“Hairdif  ■ r-A--R---  *5-4) 

"U/AIRDIF  FIT 
and 

K '“Hmorsair  ■ 12'7MOKSA1R  <5-5> 

’‘‘/ANISN  FIT 

where  the  fits  Indicated  are  least  squares  fits  of  the  AIRDIF  and  ANISN 
homogeneous  air  data  to  Murphy's  radiation  transmission  equation  (equa- 
tion (4.51)). 

Fits  were  used  to  represent  the  homogeneous  air  data  for  several 
reasons.  First,  these  fits  are  essential  to  the  direct  computation  of 
K factors  within  the  AIRDIF  and  MORSAIR  codes.  Next,  the  fits  of  the 
1-D  ANISN  data,  as  mentioned  before,  will  form  the  new  data  base  of  the 
next  version  of  the  SMAUG  code.  The  ATR  code  also  uses  these  kinds  of 
fits.  Finally,  because  of  the  statistical  nature  of  Monte  Carlo  data, 
some  fit  had  to  be  used  to  represent  the  homogeneous  air  MORSAIR  data. 

As  will  be  shown  below,  the  fits  of  the  1-D  ANISN  data  are  also  a good 
representation  of  the  homogeneous  air  MORSAIR  data. 

The  following  things  must  be  considered  In  evaluating  the  compari- 
son of  results  from  equations  (5.4)  and  (5.5): 

1.  How  well  do  the  fits  of  the  1-D  AIRDIF  and  ANISN  data  represent 
the  actual  data? 
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2.  What  type  of  statistics  are  associated  with  the  MORSAIR  data? 

3.  How  well  does  MORSAIR  agree  with  ANISN  in  homogeneous  air? 

In  general,  the  fits  to  Murphy's  transmission  equation  are  excel- 
lent representations  of  the  actual  homogeneous  air  transport  results. 

The  comparisons  presented  earlier  in  Figures  7 and  8 are  typical  of  the 
accuracy  achievable  with  these  fits.  Except  in  the  steep  build-up 
region  (0  ~ 40  g/cm2)  which  characterizes  the  secondary  gamma  4ttR2  dose/ 
range  curves,  the  fit  values  are  within  5%  of  the  actual  data.  In  the 
secondary  gamma  build-up  region,  the  fit  value?  are  generally  10-15% 
higher  than  the  actual  data. 

All  of  the  2-D  MORSAIR  data  reflect  statistical  oscillations  inher- 
ent in  Monte  Carlo  results.  For  ranges  less  than  100  g/cm2,  the  standard 
deviations  ran  between  10  and  20%  for  ranges  between  100  and  150  g/cm2 
and  between  20  and  30%  for  ranges  between  150  and  200  g/cm2.  The  scatter 
in  the  MORSAIR  K factors  is  due  mostly  to  these  statistics.  However, 
some  of  the  oscillations  clearly  lie  outside  of  these  statistics.  In 
particular,  at  mass  ranges  greater  than  about  140  g/cm2,  the  MORSAIR 
secondary  gamma  doses  seemed  to  be  biased  low. 

Figures  9-12  show  that  the  ANISN  fits  and  MORSAIR  homogeneous  air 
results  are  in  good  agreement.  In  fact,  the  primary  reason  the  ANISN 
fits  were  used  in  defining  the  MORSAIR  K factors  is  that  these  fits  also 
appear  to  be  good  fits  of  the  MORSAIR  homogeneous  air  results.  The 
biggest  differences  occur  in  the  secondary  gamma  dose.  In  the  build-up 
region,  the  1-D  MORSAIR  results  seemed  to  be  about  10%  higher  than  the 
1-D  ANISN  fits.  Also,  as  mentioned  above,  the  MORSAIR  secondary  gamma 
results  are  somewhat  lower  (10-20%)  at  ranges  greater  than  140  g/cm2. 
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Neutron  Silicon  Doses  from  a Fission  Source  in  Homogeneous  Air 


Figure  11.  Neutron  Silicon  Doses  from  a Thermonuclear  Source  in  Homogeneous  Air 
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In  summary,  all  three  of  the  factors  just  discussed  have  a hearing 
on  the  results  of  the  K factor  comparisons  and  must  be  considered  in 
analyzing  the  results.  Without  exception,  it  can  be  shown  that  any 
significant  differences  that  exist  between  the  AIRDIF  and  MORSAIR  K 
factors  are  due  to  these  factors  rather  than  differences  in  the  way  MIS 
perturbs  the  diffusion  results  and  the  Monte  Carlo  simulation  of  the 
transport  doses.  For  example,  earlier  it  was  pointed  out  that  for  ranges 
less  than  40  g/cm2,  the  AIRDIF  secondary  gamma  K factors  were  lower 
( ~ 20%)  than  the  MORSAIR  results.  This  difference  is  due  to  the  fact 
that  over  this  range  the  fit  of  the  1-D  AIRDIF  data  is  higher  than  the 
actual  AIRDIF  data,  while  the  fit  of  the  1-D  ANISN  data  is  somewhat 
lower  than  the  actual  MORSAIR  homogeneous  air  results  over  this  same 
range.  Also,  recall  that  the  AIRDIF  secondary  gamma  K factors  were 
higher  than  the  corresponding  MORSAIR  results  for  ranges  creator  than 
140  g/cm2.  This  difference  is  due  to  the  combination  of  the  2-D  air 
MORSAIR  results  being  biased  low  and  the  1-D  homogeneous  air  MORSAIR 
results  being  lower  than  the  fits  of  the  1-D  ANISN  data  over  this  range. 

In  conclusion,  comparisons  of  the  AIRDIF  and  MORSAIR  K factors  have 
been  made  using  a fission  and  thermonuclear  source,  for  silicon  and 
tissue  dose,  for  source  altitudes  between  5 and  20  km,  and  for  ranges 
out  to  200  g/cm2.  In  almost  all  cases,  the  AIRDIF  and  MORSAIR  K factors 
for  neutrons  and  secondary  gammas  are  in  good  agreement.  These  results 
support  the  validity  of  the  use  of  diffusion  theory  to  evaluate  the  mass 
integral  scaling  approximation  as  applied  to  the  atmospheric  radiation 
transport  problem.  This  result  was  not  unexpected  for  several  reasons. 
First,  it  must  be  remembered  that  diffusion  theory  is  a low  order  approx- 
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imatlon  of- transport  theory.  Comparisons  have  been  presented  in  this 
section  that  show  diffusion  theory  is  more  than  an  adequate  model  for 

t 

predicting  integral  quantities  like  dose.  Finally,  as  pointed  out  in 
Section  II,  the  mass  scaling  law  can  be  derived  from  both  the  diffusion 
and  Boltzmann  transport  equations  under  the  exact  same  set  of  assumptions. 


A Limitation 

During  the  course  of  this  research,  several  problems  were  encountered 
which  ultimately  limited  the  validity  of  the  diffusion  model  to  source 
altitudes  under  20  km  and  to  receiver  altitudes  less  than  25  km.  First, 
at  altitudes  above  25  km,  the  air  is  so  thin  (see  Figures  1 and  2)  that 
the  radiation  streams  for  great  distances  between  interactions.  The 
angular  dependence  of  the  fluence  at  these  altitudes  again  approaches 
that  of  a delta  function  which  diffusion  theory,  with  its  inherent 
assumption  of  linear  angular  variation,  is  simply  incapable  of  properly 
treating.  Use  of  the  first  collided  source  minimizes  but  does  not 
eliminate  this  problem. 

Secondly,  numerical  problems  (slow  convergence  rates,  negative 
fluences,  etc.)  were  encountered  when  AIRDIF  was  run  at  altitudes  greater 
than  25  km.  An  extensive  investigation  into  this  problem  revealed  that 
the  source  of  the  numerical  instabilities  was  the  expanding  coordinate 
system.  Specifically,  at  altitudes  above  25  km,  the  expansion  rates 
necessary  to  match  the  atmospheric  density  gradient  produce  two 
undesirable  effects.  First,  the  horizontal  mesh  becomes  too  coarse  in 
a physical  distance  sense  at  the  higher  altitudes.  This  is  particularly 
disturbing  because  although  the  mesh  interval  becomes  large  physically, 
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It  Is  still  about  equal  from  a mean  free  path  standpoint  to  the  intervals 

at  the  lower  altitudes,  where  instability  problems  were  not  encountered. 

* 

Secondly,  the  most  severe  problems  were  encountered  in  that  domain  where 
"line  of  sight"  between  the  mesh  points  and  source  point  was  interfered 
with  by  the  outer,  exponentially  curved  boundary  of  the  mesh  envelope. 
This  is  also  disturbing,  because  it  is  simply  not  possible  to  preserve 
"line  of  sight"  and  at  the  same  time  maintain  the  appropriate  mesh 
expansion  rate  at  these  altitudes  with  an  exponentially  expanding 
coordinate  system. 

Finally,  it  should  be  mentioned  that  even  if  solutions  to  the 
problems  just  discussed  could  be  found,  at  least  one  additional  problem 
associated  with  the  characteristics  of  the  U.S.  Standard  Atmosphere 
would  remain.  Specifically,  as  shown  in  Figure  2.  at  an  altitude  of 
25  km,  there  are  only  about  25  g/cm2  above.  This  means  that  no  matter 
how  high  physically  the  top  boundary  is  placed,  receivers  located  at 
altitudes  above  25  km  will  still  be  within  a couple  of  mean  free  paths 
(20-30  g/cm2)  of  the  top  boundary  where  a zero  return  current  condition 
is  applied.  The  angular  dependence  of  the  fluence  here  is  drastically 
different  from  linear.  Therefore,  the  diffusion  approach  would  have 
questionable  validity  within  a couple  of  mean  free  paths  of  the  top 
boundary  (i.e.,  above  25  km).  It  should  be  noted  that  this  problem  can- 
not be  alleviated  by  using  better  boundary  conditions.  The  zero  partial 
return  current  condition,  which  is  applied  in  AIRDIF,  is  the  best 
approximation  of  the  true  vacuum  condition  that  is  possible  to  apply 
with  diffusion  theory. 


56 


Since  the  nonorthogonal  expanding  coordinate  system  (or  one  similar 
to  it)  and  an  exponent!' al-like  atmo.-pheric  model  would  be  required  for 
any  or*der  of  numerical  solution  of  the  transport  equation  in  variable 
density  air,  it  was  felt  that  the  problems  just  discussed  probably  would 
not  be  unique  to  diffusion  theory.  Therefore,  higher  order  transport 
solutions  for  altitudes  above  25  km  were  not  pursued.  It  turns  out  that 
the  computational  efficiency  of  Monte  Carlo  techniques  improves*  at  the 
higher  altitudes,  because  only  a few  collisions  need  be  followed  before 
particles  escape  out  the  top  of  the  atmosphere.  Further,  the  Monte 
Carlo  method  is  ideally  suited  for  handling  the  streaming  radiation 
problem.  For  these  reasons,  about  six  months  prior  to  the  completion  of 
the  research  reported  here,  tht  ‘ir  Force  Weapons  Laboratory  awarded  a 
contract  to  Kaman  Sciences  Corporation  to  conduct  a high  altitude 
evaluation  of  MIS  with  the  MORSAIR  Monte  Carlo  code  (Ref.  31). 

Eamon  has  recently  completed  this  evaluation  (Ref.  10).  His 
assessment  is  based  on  the  same  sources  and  response  functions  used  in 
this  study  and  concentrates  on  source  altitudes  between  20  and  80  km. 
Eamon 's  results  show,  as  expected,  that  at  these  altitudes  leakage  is 
the  dominant  effect.  Neutron  and  secondary  gamma  K factors  as  low  as 
0.1  are  reported  for  many  cases. 


*Monte  Carlo,  even  at  high  altitudes,  is  still  too  expensive  to  run  on 
a repeated  basis  as  required  in  S/V  analyses. 
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VI . Evaluation  of  Mass  Integral  Scaling 


K Factor  versus  Mass  Range*  Plots 

The  AIRDIF  code  was  used  to  generate  silicon  and  tissue  MIS  K 
factors  for  the  fission  and  thermonuclear  sources  shown  in  Table  V.  The 
neutron  and  secondary  gamma  silicon  K factors  are  presented  as  a function 
of  mass  range  in  Appendix  E for  the  fission  source  and  in  Appendix  F for 
the  thermonuclear  source.  In  each  case,  plots  are  included  at  source 
altitudes  of  5,  10,  15,  and  20  km  and  at  seven  receiver  altitudes  for 
each  source  altitude.  These  data  were  used  to  assess  the  validity  of 
the  MIS  approximation  and  are  the  basis  of  the  discussion  below. 

Parametric  Analysis 

The  K factor  is  a function  of  mass  range,  particle  type,  source 
altitude,  receiver  altitude,  source  type,  and  dose  response  type.  In 
this  subsection,  the  sensitivity  of  the  K factor  to  each  of  these 
parameters  is  examined. 

Sensitivity  to  Mass  Range.  An  examination  of  the  data  in  the  Appen- 
dixes shows  that  over  the  5-25  km  altitude  domain,  the  K factor  is  a 
parabolic-shaped  function  of  mass  range.  For  points  close  to  the  Z-axis, 
the  K factor  is  less  than  one  and,  in  general,  is  decreasing.  Then  It 
levels  off  and  starts  to  increase  at  ranges  around  100  g/cm2.  The 
increase  continues  over  the  100-200  g/cm2  range  with  the  K factor  achieving 


*In  this  section's  figures,  the  symbol  "RH0R"  is  used  for  mass  range 
which  was  defined  in  equation  (2.3). 
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values  that  can  be  much  greater  than  one  at  large  mass  ranges. 

These  trends  are  due  to  the  combination  of  the  leakage  and  mass 
distribution  effects  described  in  Section  II.  Leakage  which  results  in 
the  scaled  doses  being  higher  is  the  primary  reason  why  the  K factors 
are  less  than  one  over  such  wide  ranges.  The  curvature  and  rise  portion 
of  the  K factor  curves  are  due  to  the  mass  distribution  effect. 

As  explained  in  Section  II,  in  variable  density  air  many  particles 
reach  receivers  at  large  mass  ranges  after  traveling  great  distances  in 
the  less  dense  air  above  the  source.  These  particles  traverse,  on  the 
average,  much  less  air  than  predicted  by  the  simple  mass  integral  and, 
therefore,  suffer  less  attenuation.  This  phenomenon  results  ir>  the  2-D 
variable  density  air  doses  being  higher  than  the  scaled  doses  at  these 
depths  and  gives  K factors  which  are  increasing  to  values  that  can  be 
much  greater  than  one.  Particle  conservation  dictates  that  if  there  is 
an  enhancement  at  these  depths,  there  must  also  be  a corresponding 
depression  in  the  2-D  air  doses  at  smaller  mass  ranges  (i.e.,  K factor 
< 1).  This  effect  is  responsible  for  the  close-in  K factors  being  less 
than  one  at  low  altitudes  where  leakage  is  insignificant. 

Sensitivity  to  Particle  Type.  Since  leakage  of  both  neutrons  and 
gammas  can  decrease  the  secondary  gamma  doses,  the  secondary  gamma  K 
factors  are,  in  general,  depressed  more  than  the  corresponding  neutron 
K factors.  This  is  illustrated  in  Figure  13  for  silicon  response  at 
coaltitude  for  a fission  source  at  20  km. 

Sensitivity  to  Source  Altitude.  As  expected,  the  depression  due  to 
leakage  increases  as  the  source  altitude  increases.  This  is  illustrated 
for  a f 1 ssi on  source  In  Figure  14  for  neutrons  and  in  Figure  15  for 
secondary  gammas. 
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Figure  14 


. Sensitivity  of  the  Coaltitude  Neutron  Silicon  K Factor 
for  a Fission  Source  to  Source  Altitude 
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SEC-GAMMA  SILICON  KF ACTOR 


Sensitivity  of  the  Coaltitude  Secondary  Gamma  Silicon 
K Factor  for  a Fission  Source  to  Source  Altitude 


Sensitivity  to  Receiver  Altitude.  In  general,  the  depression  due  to 
leakage  is  greater  for  receiver  altitudes  below  the  source  than  for  those 


above.-  Also,  the  K factor  seems  to  rise  more  sharply  for  receivers 
above  the  source.  These  effects  are  illustrated  in  Figure  16  for  neutron 
silicon  K factors  for  a fission  source  at  20  km. 

Sensitivity  to  Source  Spectrum.  As  illustrated  in  Figure  17,  at  20 
km  the  neutron  silicon  K factor  is  not  very  sensitive  to  source  type  for 
mass  ranges  out  to  160  g/cm2.  At  larger  ranges  the  fission  K.  factor 
increases  faster  than  the  thermonuclear  K factor.  The  corresponding 
comparison  of  secondary  gamma  K factors  is  made  in  Figure  18.  This 
figure  shows  that  the  secondary  gamma  K factor  is  somewhat  more  sensitive 
o source  type  than  the  neutron  K factor.  In  particular,  leakage  seems 
to  depress  the  secondary  garrma  dose  from  a fission  source  more  than  that 
from  a thermonuclear  source. 

Sensitivity  to  Dose  Response.  For  source  altitudes  less  than  20  km 
where  mass  distribution  actively  competes  with  leakage,  the  K factor  is 
not  very  sensitive  to  the  type  of  dose  response  used  in  defir, ing  the 
K factor.  This  is  illustrated  for  the  fission  source  at  20  km  in 
Figure  19  for  neutrons  and  in  Figure  20  for  secondary  gammas.  It  should 
be  mentioned  that  at  higher  altitudes  where  leakage  is  clearly  the 
dominant  effect,  Eamon  (Ref.  10)  reports  that  the  K factor  is  quite 
sensitive  to  the  type  of  response  used. 

MIS  Iso-Error  Contours 

A major  task  in  studies  of  this  type  is  to  find  a concise  and  yet 
definitive  method  of  summarizing  the  results.  The  method  developed 
here  was  to  use  the  K factor  versus  mass  range  plots  to  construct  MIS 
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Figure  T.  Sensitivity  of  the  Neutron  Silicon  K Factor  tor  a Fission 
Source  at  20  km  to  Receiver  Altitude 
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Figure  17.  Sensitivity  of  the  Coaltitude  Neutron  Silicon 
K Factor  for  a Source  at  20  km  to  Source  Type 
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Figure  18.  Sensitivity  of  the  Coaltitude  Secondary  Gamma 
Silicon  K Factor  for  a Source  at  20  km  to 
Source  Type 
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Figure  19,  Sensitivity  of  the  Coaltitude  Neutron  K Factor 
for  a Fission  Source  at  20  km  to  Dose  Response 
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Figure  20,  Sensitivity  of  the  Coaltitude  Secondary  Gamma 
K Factor  for  a Fission  Source  at  20  km  to 
Dose  Response 


iso-error  contours.  These  graphs  bound  the  MIS  associated  error  over 
various  spatial  domains  in  the  1-25  km  region.  They  were  constructed 
using  silicon  K factors  and  the  K factor/error  relationship  defined  by 
equation  (3.3).  Since  the  K factor  in  this  domain  was  found  to  be 
insensitive  to  dose  response,  they  also  apply  to  tissue  response.  The 
contours  are  defined  in  terms  of  altitude  and  horizontal  range  from  the 
z-axis. 

Contour  for  a Source  at  5 km.  The  250  g/cm2  contour  about  a source 
at  5 km  is  shown  in  Figure  21.  This  line  also  happens  to  represent  the 
10%  maximum  MTS  error  contour,  i.e.,  scaled  doses  for  receivers  inside 
Region  1 have  at  most  a 10%  MIS  error  associated  with  them.  Actually, 
most  of  the  errors  in  Region  1 are  less  than  7%.  The  errors  associated 
with  Region  2 were  not  determined  in  this  study  because  this  region 
extends  much  beyond  the  depths  to  which  the  diffusion  model  was  verified. 
This  contour  is  valid  for  both  neutrons  and  secondary  gammas  and  for 
both  the  fission  and  thermonuclear  source.  Since  the  MIS  errors  are 
smaller  at  lower  altitudes  where  the  inhomogenieties  in  the  air  are 
smaller,  the  indicated  results  are  also  valid  for  source  altitudes  between 
1 and  5 km.  At  altitudes  below  1 km,  the  air/ground  interface,  which 
was  not  considered  in  this  study,  can  significantly  perturb  the  scaled 
dose. 

Contour  for  a Source  at  10  km.  Figure  22  shows  the  220  g/cm2  and 
the  20%  maximum  MIS  error  contour  about  a source  at  10  km.  The  contour 
is  again  valid  for  both  sources,  both  responses,  and  for  neutrons  and 
secondary  gammas.  Typical  errors  in  Region  1 are  around  10%  for  neutron 
dose  and  15%  for  secondary  gamma  dose.  The  errors  associated  with  Region 
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Figure  21.  10%  Maximum  MIS  Error  Contour  for  a Source  at  5 km 
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Figure  22,  20%  Maximum  MIS  Error  Contour  for  a Source  at  10  km 
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2 were  not  determined  In  this  study  for  the  reason  cited  above. 


Contour  for  Sources  at  15  km.  At  this  altitude  the  deficiencies  of 
MIS  arfe  more  pronounced  and,  therefore,  more  details  must  be  included  in 
defining  the  error  contours.  Figure  L J shows  MIS  error  contours  for 
neutron  dose  from  a fission  source  at  15  km.  In  Region  1,  the  MIS  errors 
vary  from  0 to  25%  with  the  scaled  doses  being  higher  than  the  actual 
2-D  variable  density  air  dose.  The  errors  in  Region  2 are  also  in  the 
0-25%  range,  but  now  the  scaled  doses  are  lower  than  the  2-D  results.  In 
Region  3,  the  MIS  errors  are  in  the  25-50%  range  with  the  scaled  doses 
again  low.  Errors  greater  than  50%  are  incurred  in  Region  4.  Determina- 
tions were  not  made  in  Region  5.  The  corresponding  MIS  error  contours 
for  neutron  dose  from  a thermonuclear  source  are  shown  in  Figure  24. 

Note  that  the  fission  and  thermonuclear  source  contours  are  almost 
identical . 

The  secondary  gamma  error  contours  for  these  sources  at  15  km  are 
shown  in  Figures  25  and  26.  In  these  figures,  the  errors  in  Region  1 
va’ y from  0-33%.  In  Region  2 the  scaled  doses  are  low  with  the  errors 
ranging  from  0-25%.  In  Region  3 errors  greater  than  25%  are  incurred. 
Determinations  were  not  made  in  Region  4.  Again,  the  contours  for  the 
fission  and  thermonuclear  source  are  very  similar. 

Contours  for  Sources  at  20  km.  The  MIS  error  contour  for  sources 
at  20  km  are  shown  in  Figures  27-30.  The  fission  and  thermonuclear 
contours  indicate  that  very  significant  MIS  errors  are  easily  incurred 
at  this  altitude  because  of  leakage  and  mass  distribution  effects. 
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Secondary  Gamma  Dose  MIS  Error  Contours 
for  a Fission  Source  at  15  km 
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Figure  28.  Neutron  Dose  MIS  Error  Contours  for  a 
Thermonuclear  Source  at  20  km 
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Figure  29.  Secondary  Gamma  Dose  MIS  Error  Contours  for  a 
Fission  Source  at  2u  km 
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Figure  30,  Secondary  Gamma  Dose  MIS  Error  Contours  for  a 
Thermonuclear  Source  at  20  km 
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VII.  Summary  and  Conclusions 


The  purpose  of  the  research  reported  here  was  to  conduct  a system- 
atic evaluation  of  the  MIS  approximation  as  applied  to  the  atmospheric 
radiation  transport  problem.  The  unique  feature  of  this  research  is 
that  the  evaluation  was  conducted  with  a special  form  of  diffusion 
theory.  Ey  comparing  the  diffusion  results  to  Monte  Carlo  calculations, 
it  was  shown  that  the  diffusion  approach  is  an  adequate  method  for 
evaluating  MIS.  This  is  significant  because  the  evaluation  was  conducted 
with  diffusion  theory  at  a fraction  of  the  cost  that  would  have  been 
required  by  a comparable  Morite  Carlo  evaluation. 

The  evaluation  was  conducted  by  determining  the  K factor,  defined 
to  be  the  ratio  of  the  4trR2  dose  in  2-0  variable  density  air  to  the  ^R* 
MIS  dose,  as  a function  of  mass  range,  source  altitude,  and  receiver 
altitude.  K factors  were  determined  for  neutron  and  secondary  gamma 
doses,  for  silicon  and  tissue  response,  and  for  a fission  and  a thermonu- 
clear source.  Plots  of  the  K factor  versus  mass  range  for  source  alti- 
tudes, of  5,  10,  15,  and  km  were  prepared  and  are  included  in  the 
appendixes.  As  pointed  out  in  Section  III,  these  data  can  be  used  as  a 
multiplicative  correction  for  scaled  doses.  They  were  also  used  to  con- 
struct a series  of  MIS  iso-error  contours  which  clearly  define  the  MIS 
error  bounds  over  various  spatial  domains. 

The  results  of  the  evaluation  show  that  MIS  is  acceptable  (errors 
less  than  20%)  out  to  mass  ranges  of  220  g/cm2  for  source  altitudes 
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between  1-10  km.  At  higher  altitudes,  the  error  contour  graphs  show 
that  significant  errors  due  to  scaling  are  easily  incurred.  These 
contours  Indicate  where  an  MIS  correction  is  necessary,  and  the  K factor 
versus  mass  range  plots  can  be  used  to  extract  multiplicative  correction 
factors  for  the  scaled  doses  in  these  regions.  A significant  conclusion 
derived  from  the  evaluation  Is  that  at  altitudes  between  10  and  20  km 
the  mass  distribution  effect  completely  overwhelms  leakage  at  large  mass 
ranges.  This  results  in  scaled  doses  that  can  be  much  lower  than  the 
actual  2-D  variable  density  air  results.  Extreme  caution  should  be 
exercised  in  using  scaled  results  at  depths  greater  than  200  g/cm2. 

Finally,  a couple  of  recommendations  are  in  order.  First,  an 
accurate,  efficient,  and  validated  MIS  correction  factor  algorithm  should 
be  developed  and  incorporated  into  the  scaling  codes.  The  results  of 
a preliminary  investigation  at  the  Air  Force  Weapons  Laboratory  (AFWj  ) 
indicate  that  this  is  an  achievable  goal.  The  AFWL  study  utilized  the 
AIRDIF  K factor  data  reported  here,  and  indications  are  that  these  data 
can  be  efficiently  and  accurately  parameterized.  Secondly,  because  of 
the  importance  of  their  environments  at  high  altitudes,  it  is  recommended 
that  the  mass  integral  scaling  approximation  be  evaluated  for  prompt 
gammas*  and  X-rays. 


*As  stated  earlier,  this  study  considered  only  the  secondary  gammas 
produced  by  neutron  inelastic  scatter  with  air. 
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Appendix  A 

Derivation  of  the  Hass  Scaling  Law 

In  this  Appendix,  the  mass  scaling  law  is  derived  from  the  energy 
dependent  diffusion  equation.  The  derivation  parallels  Zerby's  deriva- 
tion of  the  law  from  the  Boltzmann  transport  equation  (Ref.  15). 

First,  consider  the  time  Independent  diffusion  equation  with  an 
isotropic  point  source  in  an  infinite  homogeneous  medium  cf  density  pA 
which  is  given  by 

Da(E)/£_  f(ra,e)  + 2— fM)  * >>>  fM 

\dRA  RA  dRA  / 


f dE*  E^(E '-+£)  f(Ra,E’)  +^A.^1sA(E) 
J S V A ' 4tR? 


3 0 


(A.  1) 


where 


Da(E)  = the  energy  dependent  diffusion  coefficient  in  units  of  cm 
F^Ra,e)  * the  energy  dependent  particle  fluence  at  a distance  of  RA 
from  the  source  in  units  of  particles/cm2-eV 
£^(E)  3 the  energy  dependent  macroscopic  total  interaction  cross 
section  at  density  pa  in  units  of  cm"1 
E^(E'-*E)  3 the  macroscopic  cross  section  for  scatter  from  energy  E' 
to  energy  E in  units  of  cm"1  per  eV 
*(ra)  = the  dirac  delta  function  (cm"1) 


S^»(£)  • the  energy  dependent  particle  source  In  wilts  of 
partlcles/eV 

Now  consider  a second  Medium,  medium  B,  composed  of  the  same  material 
as  medium  A but  having  density  pg  which  may  be  different.  In  this 
medium  the  diffusion  equation  can  be  written  as 


4(E)  f(rb,e) 


+ ‘-M  S®(E)  - 0 


(A. 2) 


The  macroscopic  cross  sections  in  medium  A can  be  expressed  as  the 
product  of  medium  B cross  sections  and  the  density  ratio  between  the 
media,  i.e., 

ZA  * Z®  ^ (A. 3) 

sj  «*  — (A.  4) 

PB 

DA  ■ Dfc  — (A. 5) 

PA 

Using  these  relationships  and  under  the  assumptions  that 

pARA  “ Pb1^  ' 

and 

SA(E)  > SB(E)  (A. 7) 

equation  (A. 2)  can  be  rewritten  as 
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, « (ra  oa^)  sA<e>  °b 

4”rA  0A 


(A. 8) 


A comparison  of  equation  (A.l)  to  (A. 8)  shays  that  they  are  Identi- 
cal except  for  the  source  term.  The  volume  integrated  source  in  medium 
A is 


I 


SA(E)  4*R*  dRA  = SA(E) 


(A. 9) 


In  medium  B,  however,  the  volume  integrated  source  is  given  by 


f ?a(e)  PB  AitR2  dR  = l[ 

sa(e) 

(A. 10) 

{ 4"ra  “a  A A V 

V 

Since  the  volume /integrated  in  medium  B differs  by  a factor  of 
from  the  source  in  medium  A,  it  must  follow  that 


or 


= (^)  ?M 

4ttR|  f(rb,e)  = 4ttRa  f(ra,e) 


(A.n) 


(A. 12) 
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where  we  have  used  equation  (A. 6)  and  included  a factor. of  4n,  as  is 
commonly  done,  to  write  the  scaling  law  in  terms  of  quantities  in  which 
geometric  attenuation  has  been  removed. 

In  summary,  the  mass  scaling  law  states  that  if  two  infinite 
homogeneous  media  composed  of  identical  material  and  containing  identi- 
cal isotropic  point  sources  are  related  by 


then. 


PB*B  “ pARA 


(A.6) 


4**1  F("B.E)  ■ 4"ra  f(rA*E) 


(A. 12) 
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Appendix  B 


The  AIRDIT  Code 


General  Description 

The  AIRDIF  code  is  based  on  the  atmospheric  diffusion  theory  pre- 
sented in  Section  IV  and  is  capable  of  computing  both  the  radiation 
environments  and  MIS  K factors  for  atmospheric  nuclear  detonations.  It 
is  a structured  computer  code  composed  of  42  independent  subroutine 
modules  and  written  in  ANSI  standard  FORTRAN. 

Figure  B-l  shows  the  major  hierarchy  of  the  AIRDIF  code.  The  logic 
in  a structured  code,  like  AIRDIF,  always  flows  from  the  top  downward. 

For  example,  the  executive  routine,  EXEC,  calls  three  routines:  DATAIN, 

CALC,  and  RESGUT.  EXEC  first  calls  DATAIN  to  read  in  the  input  data. 
Control  then  shifts  back  to  EXEC  where  subroutine  CALC  which  performs 
all  of  the  required  calculations  is  called.  Control  again  shifts  back 
to  EXEC,  and  RESOUT  which  outputs  and  stores  the  results  is  called. 

There  are  many  advantages  inherent  to  the  characteristic  modularity 
of  such  structured  programs.  While  it  is  beyond  the  scope  of  this 
appendix  to  enumerate  on  all  of  them,  one  in  particular  is  worthy  of 
mention.  Specifically,  such  codes  are  inherently  very  versatile.  For 
example,  the  onl;  routine  which  is  unique  to  diffusion  theory  is  MATRIX 
which  defines  the  finite  differenced  diffusion  equations.  A higher  order 
method  such  as  Roberd's  synthesis  method  (Ref.  32)  could  be  incorporated 
into  AIRDIF  by  simply  replacing  this  routine  with  one  based  on  the  finite 
differencing  of  the  synthesis  equation. 


Figure  B-1.  Hierarchy  of  the  AIRDIF  Code 


AIRDJF  has  been  successfully  run  on  both  the  CDC-6600  and  7600 


■uto*-.. 


computers.  Execution  times  of  about  0.05  sec  per  mesh  point  (for  37  neutron 
and  2V  gamma  groups)  are  required  on  the  7600.  Since  most  problems  require 
somewhere  between  1500  and  3000  mesh  points,  a couple  of  minutes  of  7600 
CP  time  is  normally  required  for  each  run.  In  its  current  form,  AIRDIF 
has  storage  set  aside  for  up  to  5000  mesh  points  and  requires  approxi- 
mately 33  octal  K of  small  core  memory  (SCM)  and  220  octal  K of  large 
core  memory  (LCM)  on  the  7600. 

Computation  of  Radiation  Environments 

AIRDIF  utilizes  the  analytical  first  collided  source  defined  by 
equation  (4.31)  and  numerically  solves  equation  (4.50)  group  by  group 
to  obtain  the  collided  fluence  at  all  mesh  points.  Total  fluences  in 
each  group  are  then  obtained  by  combining  the  uncollided  and  collided 
fluences.  User-supplied  dose  response  functions  and  equation  (3.1)  are 
then  used  to  compute  the  dose  at  all  mesh  points. 

As  pointed  out  in  Section  IV,  AIRDIF  is  capable  of  determining  the 
radiation  environments  in  both  homogeneous  air  and  in  the  U.S.  Standard 
Atmosphere  (2-D  air).  The  user  simply  specifies  the  option  desired. 

For  homogeneous  air,  the  conditions  specified  by  equations  (4.36) 
through  (4.38)  are  applied  internally  and  the  resulting  five-point 
difference  equations  are  solved  iteratively.  In  2-D  air,  a nine-point 
difference  equation  is  solved. 

There  is  also  an  option  in  AIRDIF  for  storing  all  results  including 
the  differential  energy  fluence  at  each  mesh  point  on  a permanent  file. 

This  data  file  can  then  be  accessed  later,  and  environments  in  any  dose 
units  can  be  quickly  computed  with  a simple  editing  code. 
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Computation  of  K Factors 

AIRDIF  is  also  capable  of  directly  computing  K factors  defined  by 
equation  (3.2).  An  AIRDIF  homogeneous  air  run  is  made  first,  and  the 
results  for  a particular  source  are  stored  on  a permanent  file.  An  edit 
code  called  MISFIT  (Ref.  29)  is  then  used  to  access  this  file  and  to  per- 
form a least  squares  fit  of  the  dose  for  a user-supplied  dose  response 
function  to  Murphy's  radiation  transmission  equation  (Ref.  16:36), 

T(x)  ■ exp^A  + Bx  + Cx2  +Dx3/2  + Ex1//a  + Fx1/3  + G£n(x)j  (4.51) 

where 

T(x)  = 4uR21-D  MIS  dose  at  a mass  range  of  x g/cm2. 

The  fit  coefficients  for  both  neutrons  and  gammas  are  included  as 
part  of  the  input  for  the  2-D  variable  density  air  run.  Using  these 
coefficients  and  equation  (4V51),  AIRDIF  computes  the  4irR2  1-D  MIr- 
doses.  These  results,  the  calculated  2-D  air  doses,  and  the  K factor 
equation  (3.2)  are  then  used  by  AIRDIF  to  compute  the  K factors  for  the 
specified  source  and  dose  response  at  all  mesh  points. 


Appendix  C 


Radiation  Dose  Comparisons 


This  appendix  presents  graphs  which  conpare  the  diffusion  radiation 
environments  to  transport  theory  and  Monte  Carlo  environments.*  In 
homogeneous  air,  the  AIRDIF  diffusion  results  are  compared  to  ANISN 
discrete  ordinates  results.  In  2-D  air  (U.S.  Standard  Atmosphere),  the 
comparisons  are  made  against  MORSAIR  Monte  Carlo  results.  The  compari- 
sons include  results  for  neutron  and  secondary  gamma  environments,  for 
silicon  and  tissue  dose  response  (Tables  III  and  IV),  for  a fission  and 
thermonuclear  source  (Table  V),  and  out  to  mass  ranges  of  200  g/cm2. 


♦Note  that  4ttR2  dose  is  plotted  as  a function  of  mass  range  in  all  graphs. 
Mass  range,  <pR>,  is  defined  by  equation  (2.3).  In  this  appendix,  the 
symbol  "RHOR"  is  used  for  <pR>. 


94 


t 


Figure  C-2.  Secondary  Ganna  Silicon  Doses  from  a 
Fission  Source  In  Homogeneous  Air 
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Figure  C-4.  Secondary  Gamma  Silicon  Doses  from  a Fission  i 

Source  at  15  km  and  Receivers  at  11.3  km  ] 
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Figure  C-6.  Secondary  Gamma  Silicon  Doses  from  a Fission 
Source  at  15  km  and  Receivers  at  15.29  km 
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Figure  C-7,  Neutron  Silicon  Doses  from  a Fission  Source 
at  15  km  and' Receivers  at  21.66  km 
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Figure  C-8.  Secondary  Gamma  Silicon  Doses  from  a Fission 
Source  at  15  km  and  Receivers  at  21.66  km 
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Figure  C-9 
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Neutron  Silicon  Doses  from  a Fission  Source 
at  20  km  and  Receivers  at  20.61  km 
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Figure  C-12,  Secondary  Gamma  Tissue  Doses  from  a Fission 
Source  at  20  km  and  Receivers  at  20.61  km 
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Figure  C-13.  Neutron  Silicon  Doses  from  a Thermonuclear 
Source  at  20  km  and  Receivers  at  20.61  km 
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Figure  C-14.  Secondary  Gamma  Silicon  Doses  from  a 
Thermonuclear  Source  at  20  km  and 
Receivers  at  20.61  km 
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Appendix  D 

i 

K Factor  Comparisons 


This  appendix  contains  comparisons  of  the  AIRDIF  and  MORSAIR  K 
factors.  These  comparisons  are  listed  in  Table  VI  and  discussed  in 
Section  V.  K factor  comparisons  are  made  for  neutrons  and  secondary 
gammas,  in  terms  of  silicon  and  tissue  dose  response,  for  a fission  and 
thermonuclear  source,  and  out  to  mass  ranges*  of  200  g/cm2. 


*Mass  range,  <pR>,  is  defined  by  equation  (2.3).  In  t;-:  appendix,  the 

symbol  "RHOR"  is  used  for  <pR>. 
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Figure  D-l.  Neutron  Silicon  K Factors  for  a Fission 
Source  at  5 km  and  Receivers  at  5.07  km 
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Figure  D-2.  Secondary  Gamma  Silicon  K Factors  for  a 
Fission  Source  at  5 km  and  Receivers  at 
5.07  km 
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Figure  D-4. 

Secondary  Gamma  Silicon  K Factors  for  a 
Fission  Source  at  10  km  and  Receivers 
at  8.84  km 
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Figure  D-5.  Neutron  Silicon  K Factors  for  a Fission 
Source  at  10  km  and  Receivers  at  10  km 


Figure  D-6,  Secondary  Gamma  Silicon  k Factors  for  a 
Fission  Source  at  10  km  and  Receivers 
at  10  km 
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Figure  D-8.  Secondary  Gamma  Silicon  K Factors  for  a 
Fission  Source  at  10  km  and  Receivers 
at  1?  13  km 
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Figure  D-ll.  Neutron  Silicon  K Factors  for  a Fission 
Source  at  15  km  and  Receivers  at  15.29  km 
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Figure  0-12.  Secondary  Gamna  Silicon  K Factors  for  a 
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Figure  D-13.  Neutron  Silicci  K Factors  for  a Fission 
Source  at  15  km  and  Receivers  at  21.66  km 


Figure  D-14.  Secondary  Gamma  Silicon  K Factors  for  a 
Fission  Source  at  15  km  and  Receivers 
at  21.66  km 
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Figure  D-17.  Neutron  Silicon  K Factors  for  a Fission 
Source  at  20  km  and  Receivers  at  20.61  km 
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Figure  D-19.  Neutron  Silicon  K Factors  for  a Fission 
Source  at  20  km  and  Receivers  at  22.08  km 


Figure  D-20.  Secondary  Gamma'  Silicon  K Factors  for  a 
• Fission  Source  at  20  km  and  Receivers 
at  22,08  km 
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Figure  D-21.  Neutron  Tissue  K Factors  for  a Fission 

Source  at  20  km  and  Receivers  at  14.85  km 
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Figure  D-22.  Secondary  Gamma  Tissue  K Factors  for  a 
Fission  Source  at  20  km  and  Receivers 
at  14.85  km 
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Figure  D-23.  Neutron  Tissue  K Factors  for  a Fission 

Source  at  20  km  and  Receivers  at  20.61  km 
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Figure  D-24.  Secondary  Gamma  Tissue  K Factors  for  a 
Fission  Source  at  20  km  and  Receivers 
at  20,61  km 
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Figure  D-27.  Neutron  Silicon  K Factors  for  a Thermonuclear 
Source  at  20  km  and  Receivers  at  14.85  km 


Figure  D-28.  Secondary  Gamma  Silicon  K Factors  for  a 

Thermonuclear  Source  at  20  km  and  Receivers 
at  14.85  km 


123 


SEC-6MBU  SILICON  KFACTCR  KE’JTRON  SILICON  KF ACTOR 


4.0 


3.0 


2.0 


1.0 


0.0 


— i — i — [ — i — i — i — i — | i — i — i i — i — i — n r 


D 


OAIRDIF 

□ MORS A I R O 

o 

O o 


O °Q 


u o 


fe><£)  cP^o  ^cPq  Qo  8 0 ^ a a 


J I L 


J I I 1 1 I I I 1 1 1 ! 1 1 ! — J- 


50 


ICO 


RH0R  (G/cir.2) 


150 


200 


Figure  D-29.  Neutron  Silicon  K Factors  for  a Thermonuclear 
Source  at  20  km  and  Receivers  at  20.61  km 
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Figure  D-30.  Secondary  Gamma  Silicon  K Factors  for  a 

Thermonuclear  Source  at  20  km  and  Receivers 
at  20.61  km 


— | — | — t — , — | — , — , — i — i — j — i — i — i r~ 

"T  * 1 

i — r 

•M 

* 

o - 

OAIRDIF 

o 

^ a M0RSAIR 

o 

— 

- 

o 

- 

- 

o 

. 

o 

— 

— o 

— 

- o o 0 

"°o  0 o^cP'a  °cp  O 861  ® °a°o  o ° 

a 

D 

o- 

L_J 1 L_J 1 1 1 1 1 1 1 L_ 

1 1 !_ 

J — l — 

ji  SO  100  11)0  200 

RH0R  (G/cir,2) 


iaiagiaMdailai 


124 


|"1V  J^nn. 


*rr<  pw-yyt-.,  VTS.,T, 


V 


j 


i 


Appendix  E 

K Factor  versus  Mass  Range  Plots  for  a Fission  Source 


This  appendix  presents  a series  of  K factor  versus  mass  range* 
plots  for  the  fission  source  in  Table  V.  The  K factors  are  based  on  the 
silicon  dose  response  functions  shown  in  Tables  III  and  IV.  Plots  are 
included  for  neutrons  and  secondary  gammas,  at  source  altitudes  of  5, 

10,  15,  and  20  km,  and  at  seven  receiver  altitudes  for  each  source  alti- 
tude . 


Index 

Run  No.  Source  Altitude  (km)  Page 

1 b 127 

2 10  134 

3 15  141 

4 20  148 


*Mass  range,  <pR>,  is  defined  by  equation  (2.3)  in  Section  II.  In  this 
appendix,  the  symbol  "RHOR"  is  used  for  <pR>. 
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Appendix  F 


K Factor  versus  Mass  Range  Plots  for  a Thermonuclear  Source 


This  appendix  presents  a series  of  1C  factor  versus  mass  range* 
plots  for  the  thermonuclear  source  in  Table  V.  The  K factors  are  based 
on  the  silicon  dose  response  functions  shown  in  Tables  III  and  IV. 

Plots  are  included  for  neutrons  and  secondary  gammas,  a*  source  alti- 
tudes of  5,  10,  15,  and  20  km,  and  at  seven  receiver  altitudes  .or  each 
source  altitude. 
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*Mass  range,  <pR>,  is  defined  by  equation  (2.3)  in  Section  II.  In  this 
appendix,  the  symbol  "RHOR"  is  used  for  <pR>. 
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